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NUMERICAL STUDY OF THE STABILITY OF A HEATED, 
WATER BOUNDARY LAYER 


Abstract 
by 


ROBERT L. LOWELL, JR. 


Numerical sclutions are obtained to the sixth order system of 
disturbance equations describing the hydrodynamic stability of a 
laminar, heated, flat plate water boundary layer, including all mean 
and disturbance property variations. The results are compared with 
those calculated using the fourth order system of Wazzan, Okamura, 
and Smith [4] (which assumes only mean viscosity variation). 

Bota systems predict significant boundary laver stabilization 
‘increased minimum critical Reynolds Number, decreased disturbance 
amplification rates, etc.) with moderate heating, but display a 
maximum and subsequent decrease as the plate surface-to-free stream 
temperature difference is further increased. Over the normal liquid 
range of water, the results of the sixth order calculations show only 
a slight enhancement of stability over those obtained from the fourth 
order system, Apparent gains in stability predicted when only dis- 
turbance viscosity terms are included in the calculations are offset 
by the further inclusion of density fluctuation terms. 

Tne insensitivity of stability characteristics to the addi- 


Plone! esnelderation of thermal disturbances is attributed to the 
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very limited region in which such disturbances are important for 


fluids of large Prandtl number such as water. 
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CDAPTER I 


INTRODUCTION 


To gain some insight into the mechanism of laminar-turbulent 
boundary layer transition, considerable experimental and theoretical 
effort has been expended in studying the hydrodynamic stability of 
a laminar boundary layer. This work has revealed that alterations 
to the mean flow by pressure gradients, boundary roughness and 
compliancy, suction or blowing, heating and the like can cause 
significant changes in stability characteristics of the flow. The 
present study focuses on the influence of heating on the stability 
of a liquid boundary layer, and is motivated by the favorable effect 
that surface heating may have on drag reduction of submerged under- 
water vehicles, 

An indication of the qualitative effect of surface heating 
can be obtained by utilizing the inviscid point-of-inflection 
criterion (for constant density flows velocity profiles having a 
point of inflection are unstable). At the surface of a flat plate, 
the boundary layer equations for a fluid flow with a nonconstant 


viscosity reduce to 
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If the plate is heated (such that TI» then (dT/dy), «0. Because 
the viscosity of a liquid decreases with increasing temperature 
(du/dT « 0), then (éu/dy).. » 0. Conversely, for a cooled plate 


Eu < To)» (du/dy) < 0. As both u, and (du/dy) are positive for 


o> 


y 
MO 


either case, it follows directiy that for T , Sr 


a is vanishingly small but negative 


Since the curvature, ay 
as y + ©, it follows that if its value at the wall is positive, there 
must be at least one point of inflexion in the boundary layer (i.e., 
some point at which the curvature is zero). Thus it is anticipated 
that heating should stabilize and cooling destabilize the liquid 
laminar boundary layer. This conjecture has been verified 
theoretically by a number of authors. (Note that it is just opposite 
to the conclusion reached for gases where viscosity increases with 
temperature). 

DiPrima and Duan [1] cite unpublished numerical results of 
McIntosh which indicate that for the two-dimensional Tollmien- 
Schlichting type of instability in a heated water boundary layer, 


the minimum critical Reynolds number, Re is increased 


Mins crilt.” 
. O e 
tenfold over the isothermal case for a 50 C temperature diiference 
o 
EE ER, 
CO 
Hauptmann [2] utilizes iutepral methods and a perturbation 
procedure to correlate variation in a velocity profile shape factor, 
arising from small changes in a temperature dependent viscosity, 


with the stability characteristics for a class of single parameter 


velocity profiles. llis results also predict very strong stabilization 





La? 


in water for relatively small heating, and are quite compatible 
with those of McIntosh cited above. 

The most extensive study of the stability of a heated and 
sed water boundary layer with a Falkner-Skan flow has been con- 
ducted by Wazzan, Okamura, and Smith [3,4 E and Wazzan, Keltner, 
Okamura, and Smith [6]. For the disturbance equation, assuming 
explicit fluid property variation in only the mean viscosity and 
no fluctuating temperature field, they obtain for all their 


analyses a modified Orr-Sommerfeld equation 
E zm a ah 
(AB A) UD 
E A gie stéi + a (0-0) 0a 
(oo) 


"here denotes differentiation with respect to y/6, In con- 


where 


trast, the variable mean-flow coefficients are obtained by solving 
the coupled mean momentum and energy boundary layer equations in 
which all fluid properties are assumed to vary [7]. In general, 
their results show that while cooling the wall does destabilize 

F ’ E | 
the flow, moderate heating very strongly stabilizes it. However, 


as heating is increased, Re ., reaches a maximum with T and 
min,crit, y 


then decreases. They explain this bchavior as follows [5]: 


e ER ee ee ae ee ee a rn te ee 


While the resuits of reference [3] are qu2litatively correct, 
they are quantitatively inaccurate duc to an error in the numerical 
intesration of equation (1.2) and subsequently have been corrected 


1,5]. The first reference still provides the best explanation of 
the numerical procedure used, 





"AS T, is increased above T. , the velocity profile becomes 
more stable, whereas the variable viscosity term, u ( n), in the 
modified O-S equation tends to destabilize the flow. However, at 
moderated rates of heating, the effects of the mean velocity pro- 


files u (n) ard v" (n) upon Re are the dominant ones. As 


Min. crite. 


the rate of heating increases, the negative effect of u (n) upon 


Ke ; increases whereas the effect of heating upon u (n) and 
HA. CELTA. i 


u" (n) and, in turn, their respective stabilizing effects on 


begin to level off, Re. h then reaches a maximum 
mi crit: 


` 


mun.crit. 


at (T) .. where the effects of y (n) and u (m) and u" (n) upon 
7 cr]L, 
ec. f balance out. As T, is increased beyond (T)  ,, , the 
Eum GS. VW v CLIL: 


destabilizing effect of u (n), becomes dominant and EN E 


decreases upon increasing T, beyond D : 


Similar results and conclusions are obtained by Potter and 
Graber [8] for the stability of Plane Poiseuille water flow with 
heat transfer. Using for their analysis the same equation as 
Wazzen, et.al. and neglecting all property variation other than 
viscosity in the solution of the mean equations, they found that 
omission of the mean viscosity gradient terms from equation (1.2) 
jmdieated inereased stability with heating. In contrast, inclusion 
of these terms show destabilization with increased heating (i.e., 
a decrease of Boun, ein Ji 
The motivation for the present investigation stems from an 


attempt to more fully explain the stability anomalies found by 


Wazzan, Okamura, and Smith through a more complete reformulation 





of the problem. Specifically, it is unclear why those authors 
solve the coupled mean flow system of momentum and energy equations 
to provide the variable coefficients for the disturbance equation 
(1.2), and yet neglect a disturbance energy equation coupled with 
the momentum equation through fluctuating fluid property variation. 
As an a priori estimate of the relative magnitude of all disturbance 
quantities and their derivatives is impossible throughout the 
entire boundary layer, the omission of temperature fluctuations in 
previously cited investigations is seemingly arbitrary and without 
justification in a situation involving heat transfer. Accordingly, 
the complete linearized, parallel flow disturbance equations (ex- 
cluding only negligible expansion work and viscous dissipation) 
eventually will be formulated and solved. 

A simple illustration substantiates the preceeding conten- 
tion. Assuming for the moment that the only property variation of 
importance is that of viscosity, it will be shown that the stability 


equations can be written more completely as 
; A Ü m Hil 2 7j !! E 
(g-cY d'-« d) -óU = -æa (9 "Zi ro a 


+22 (0"- 8°) +a" (g « o'£)] - 





— 
edi ¡(may + oec m | 
ot Ke 





coupled vith a disturbance energy equation 


SE : T mm Ir: (1.4) 
en c)+H9 ait ed 


through the relation between the viscosity fluctuations, m, aná 
the temperature fluctuations t. Assuming further that disturbance 
quantities, mand €, as well as their derivatives, m" and $6", are 
each of comparable magnitude at some point in the boundary layer 
and realizing JUS is of the same magnitude as yu", then terms 
arising due to the presence of the disturbance viscosity such as 

wo N 
mu and e mu are likewise of comparable magnitude to terms 
retained by Wazzan, et.al.such as div" ALTE , and so should 
not be discounted in the analysis. Since Wazzan, et.al. attribute 
much of the stabilization to the mean viscosity gradient terms eile 
it would seem as if omission of these terms is "not well justified" 

^J 

and it is not permissible categorically "to assume that T = 0 
throughout the boundary layer" [3]. 

Thus the focus of the investigation is upon establishing the 
difference between the fourth order system of Wazzan, et.al., 
equation (1.2), and a more complete version of the sixth order 
system represented by equations (1.3) and (1.4) in which all 
property variation, mcan and disturbance (acceptable within the 


locally parallel flow and boundary layer assumptions), is 


considered, 





Sch 


A complete discussion of the assumptions made in formulating 
the problem considered in this investigation is given in Chapter 
Il. The assumed fluid property-temperature relationships, numerical 
techniques, and assessment of the numerical errors incurred in the 
process are discussed in Chapter III. The results of this analysis 
are presented in Chapter IV, and finally the conclusions drawn from 
these results are discussed in Chapter V. The Appendices are in- 
cluded to present relevant but nonessential material that would 
otherwise obscure the primary thrust of this investigation in the 


text, 


A 
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CHAPTER II 


FORMULATION OF THE STADILITY PROBLEM 


This investigation examines the spatial stability to small- 
scale disturbances of the steady, laminar, boundary layer flow of 
a viscous, heat-conducting liquid. Since formulation of the 
problem can easily be made with far greater generality than is 
eventually required, such a procedure is foljowed to demonstrate 
that within the assumptions to be made, the additional complexities 
fail to significantly alter the resulting stability equations. 
Thus while the problem is posed for general Falkner-Skan flows 
Cu, = SW with constant cross flows (wr. = constant) and arbitrary 
disturbance orientation, solutions discussed in Chapters VI and V 
are specifically for two-dimensional, flat plate boundary layers 
with two-dimensional disturbances. The Falkner-Skan flovs are 


often referred to as "wedge flows" since for M > O they correspond 


to a uniform flow over a yawed, rigid, infinite wedge with included 


dach 
angle E [9, p 141] at zero angle of attack in the neighbor- 


A 


hood of the stagnation point. 


2.1 Basic Assumptions and Governing Equations 
To describe this flow field mathcmatically, certain basic 


assumptions are made concerning the nature of the fluid involved. 


(1) The characteristic scales for mean motion and 





H 


disturbances to 1t must be significaatly greater than those 
associated with molecular motions in the fluid (continuum 
hypothesis). This permits ascription of meaningful values to 
‘fluid and flow parameters "at a point". Further, for disturbance 
periods (2x/w) much greater than the molecular time scale, say 
the fiuid's structural relaxation time (T .) Cie. wT << JOE 
flow can be considered to be in statistical equilibrium so that 
molecular transport effects can be represented by transport co- 
ellicients Gi, k, otc.). 

(2) Once the continuum hypothesis is established, to define 
the stress relations between fluid elements, it is assumed that 
the fluid is Newtonian (linear stress variation with rate of strain) 
and isotropic (no intrinsically preferred elemental spatial 
orientation). 

These two groups of assumptions form the foundation for the 
Navier-Stokes equations (2.1-2.3). Since there is no experimental 
evidence to suggest their inadequacy under ordinary conditions, 
they will be used without reservation in this analysis, 

(3) Assuming the fluid to be a liquid, the form for the 
state equations is defined. The relative pressure independence of 
all fluid properties permits their expression explicitly as 


functions of temperature alone. 





10 


Thus, (dimensionally) 


Continuity: Bg > Qe. m (2:15 
Momentum: DUE $ SE OTe Ce 
CDE L ; dE 


Ene : DT ; dp (253 
p eot! Er + AE +2 us + SS ) 


^r ^R 


State: ES Co, <A rn = Maa ee (2.4) 
D ò d 
where — = material derivative = PA 
Dt eb Jd, 
J 
@i = ect reserernsan = =o o 
y n | SC | KK oO 
E pra te=ot=s train tensor = St OU: a Qu. d: 
oF. SE 
Í 
9, = body force 
= volumetric expansion coefficient = (E ) 
C S 
PUR = "mechanical" - thermodynamic pressure 
difference = - ee £ NG 
~ KK 
This last equation relating the mechanical pressure (negative 
of the mean normal stress, - 0,0) and thermodynamic pressure is 
a INN 


valid because only the near past is considered (or wT << yeas 
previously discussed [10] ). Note that cquations (2.2) reduce to 


8 » s À 2 e 
the Navier-Stokes equations if E. -— —y (Stokes! hypothesis, 
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valid only for an ideal monatomic gas) or Ey = 2 
For Liquid in a heated boundary layer, there is no justifi- 
cation for assuming either statement is valid in general. 


am 
EIS 
u 


‘Liebermann [11] has shown that for a number of liquids, 
positive, considerably greater than unity in some cases, and wideiy 
different for different liquids. Furthermore, for low frequencies 
CM «« 1), Narasimham (12) notes that the viscosity ratio is 
independent of disturbenee frequency. and almost independent of 
temperature [13]. (For water, Pinkerton [14! has shown a — 
variation of about 16% over the range 0-60?C as compared with a 

y variation of 118% over the same range - see Figures B.2 and B.6). 
Alternatively, for flows with heat transfer (particularly cases 
with large thermal gradients), when the density is identified as a 
function of temperature, it will not be constant in that flow 


field. In such a situation, then Sek will vanish only at a rigid 


surface for steady flovs. Thus if relatively large amplitude 





density fluctuations were expected in a liquid for which was 
large (e.g., for benzene, this ratio is 0 (10?) [11], then there 
could be a large disparity between mechanical and thermodynamic 
pressures, It will be shown in Appendix A and Section 2.3.1 that 
for values of E: order unity, the stability problem can be 
solved independently of this quantity. 

When the flow is characterized by a single time scale, namely 


that associated with the frequency of the small, rapid perturbations 


on the basic flow, each of the flow and transport parameters can 





n 


be constructed as the sum of a time-averaged mean and a fluctuating 


component, 


Qi pus Qu (Es) S Qt) (2.6) 


Equations describing the behavior of these components are developed 


in sections 2.2 and 2.3, respectively. 


2.2 The Mean Flow 

The mean equations for boundary layer flow over a wedge can 
be deduced from equations (2.1-2.3) (see Appendix A) after choosiug 
the orthogonal coordinate system shovn in figure 2.1 and making the 
following simplifying assumptions: 

(1) The basic flow is steady but not necessarily spatially 
homogeneous (eqn. 2.6), so the mean eauations evolve from a time- 
everage of the instantaneous equations. 

(2) In the momentum equations for motion parallel to the 
surface, viscous and inertial forces are of the same order of 
magnitude. Using the scaiing resulting from this requirement and 


l - . 
) or smaller, boundary layer equations 


eliminating terms of 0 (Res 
are obtained, 

(3) The wedge is "infinite" in extent in the z direction, so 
all mean flow variation is independent of the z coordinate. 


(4) Gravity is the only body force present, acting in a 


plane normal to the z axis at an angle xy to the y axis. 





< 





Coordinate system for the stability 


Fipure 2.1 
layer flow. 


4 
* 


direction of wave 
» propagation 


sce reference [15] 


analysis of a general three-dimensional boundary 


E 





Additional simplifications leading to exclusion of specific 
effects on the flow (e.g., dissipation, compression energy, buoyancy) 


will be made after examining relevant nondimensional coefficients. 


L2 Mean Flow Equations 
Using the assumptions discussed in section 2.2, equations 


(2.1-2,3) reduce to the (time-averaged) mean boundary layer 


equations (dimensionally) 


Continuity: ou EE O (25) 
AX y 
Momentum: 5 (à ou +U E A ed sk - dpi (2.8) 
( AN y : EK 
ò 4 95) 
"as, 


OB e (2.9) 


ar) 


Energy: us (a oT +V gc SI LI do (2.11) 


Before solving, several modifications are made to the x 
momentum equation. Evaluation of this equation at the boundary 
layer edge reveals the nature of the mean pressure variation im- 


pressed on the boundary layer (dimensionally) 
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pm to a, E Qe d Sin X (2.12) 


Further, due to the state equation (2.4), the density can be ex- 


panded in a Taylor series in temperature, so dimensionally 


Zeit lo + = RES E PAS (de ). m 


mcs 2 T 


or nondimensionally, using the volumetric expansion coefficient, 


B, defined in equations (2.5), 


lcg CD Qu Too (T- Cs 





ISS 


(such as is the case for 


co 
water with moderate temperature differences), this can well be 


Lor small values oí ES jo and Tæ- 1 


approximated by 


| - e Cr) = Qu To (T- 1) (2.13) 


Thus, with the modifications indicated by equations (2.12) 


and (2.13), the x momentum equation finally becomes (dimensionally) 





e mm 


^70 
fos 
E 
C^ KV 
pue 
s 
<i 
oY 
cl 
DN 


LL d 


e E Ss 
Goo Ue SES + Ro Rog SiINX (T p 


n d Lo T (2.82) 
OY OY 
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Since all fluid properties are assumed to be variable, 
particularly density, a modified Howarth-Dorotnitsyn transformation 
is suggested [7], reducing the boundary layer equations to a 
nearly incompressible form. The required independent variables 


are 


q fa e dy C2714) 


Satisfaction of the continuity equation is accomplished by intro- 


T. 


Ze 
ducing a dimensional stream function p such that 
* 
A 
psc 
D 


=*=% DE 
V nn 
S Ox 





and a nondimensjonal stream function F so that 


x » l (2.15) 


E 
pl 


where the two are related by 


(Rohe uk) E 


Defining the displacement thickness 


Pu 20) dy (2.16) 
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or in transformed variables 


E M E Cries 
vet Hd 
Ó 
hen VE ` S (2.17) 


where 


2 In - n en , (2,18) 


a constant for specified wall and free stream temperature and wedge 


angles. 


Additionally, defining 








DE o + Te 
W* = We LG Col 
then equations (2.7,2.8a,2. 10, and 2.11) transform to 
"I UM m (EYEE Een c 
i | Re E 
(ac) Gi ar (=) 
E E m m 0 
GET) + | (Mg) + RE iga] 
Je S pe KI Ge v 
+ (=) E ]* i» = Sin X (2.19) 
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EA 
= Ó 
where Re = Reynolds Number > Mc 
Es 3 
e $ 5 I 
Gr = Grashof Number = Oo 9 ar SE Ô 
We ia 
Ue EN 
E = Eckert Number = PN 
Peo | eo 
Pr = Prandtl Number = Cpe Meo 
Se 
M = E due = constant since u ids of the 
Ue dx e 


form in E for the wedge flov. 


At this point an assessment is made-of RS various force and 
energy terms cf the above equations to determine their influence 
in establishing the mean velocity and temperature boundary layer 
profiles, with the following conclusions. 
(1) Buoyancy effects can be neglected if = sin A 
This is easily accomplished analytically by requiring 
(a) the surface normal be parallel with the gravity 
vector so that y = 0, or 
(b for large x angles (i.e. sin y = 00), Ua be 


—À 


AE Gr : 
sufficiently large that pocas ieee tec ens 
: 2 
parametric group varies inversely with un: Thus, 
except for surface orientations nearly parallel to 
the gravity vector and low speeds, free convection 
can be neglected. 
(2) Dissipation and expansion (or pressure) work are 
negligible, requiring that Pr E << 1. (An estimate of these terms 


for a water flow indicate that they can never be significant while 


the boundary layer approximation remains valid). 
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As a consequence of these assumptions, the boundary layer 
equations (2.19) are explicitly independent of x. The original 
set of partial differential equations has been transformed to one 
of ordinary differential equations, so the velocity and temperature 
profiles are given by similarity solutions. 


The mean flow equations can be written finally as 


e ED + M (4 - F^) + (MER O 








New 9 (2.20) 


Note that the z momentum equation is not coupled with the 
other two; that is, F and H are determined independently of G 


provided DUE and k are also independent of G, 


2.2.2 Mean Flow Boundary Conditions 


The boundary layer equations in section 2.2.1 are 
supplemented with the following boundary conditions: 
aX =m 
at y = 0 u= y=0 (no slip at the wall) 


=x% 2 
v e O (impermeable wall surface) 


T = T (constant wall temperature) (2.21) 





20 


—X y 
ag y > œ u EL 


= 
war" 
e 


N 
od 
co 


Note that although u, may vary with x due to a nonzero 
wedge angle, the free stream temperature (identical with that at 
the boundary layer edge) is assumed constant. 

In terms of the transformed variables, these conditions 


become 


at n = 0 P nz CG = 0 
F = 0 
GE (2.212) 
as nz o F* > 1 
G >1 
H >0 


2.3 The Disturbancé Flow 

The disturbance equations are obtained as the difference be- 
tween the instantaneous and mean equations (see Appendix A). 
Simplification is accomplished by the following considerations. 

(1) Consistent with the assumptions of section 2.2, 
“infinitesimal” disturbance amplitudes permit linearization of the 
disturbance equations. (Note that all terms of this order vanished 
in taking the mean). This of course eliminates terms with 


correlated disturbance quantities (as for the mean equations) and 
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reduces the disturbance equations to a set of linear, partial 
differential equations with nonconstant TOEN flow coefficients. 
(2) The mean boundary layer flow is assumed to be locally 
‘parallel, so that for a specified distance from the leading edge, 
x, these coefficients of the disturbance quantities are functions 
only of the distance normal to the solid surface, y. That is, 
the stability of a velocity profile at a given x is to be con- 
Bee independently of the rest of the boundary layer. The 
requirements necessary to make such an assumption are elucidated 
in Appendix A and briefly summarized below. When the amplitude of 
all disturbances and derivatives of disturbances are equal (though 


not necessarily equal to each other), comparison of the mean flow 


—$ —* 
coefficients in a boundary layer sense indicates that v << u , 
— 3Q 2Q E 
w and << . Furthermore, if disturbance wave lengths 








ox oy 


are small compared with lengths characterizing changes with x of 
mean quantities, then these quantities can be considered to be 
slowly varying functions of x compared with the fluctuations. 
Again, the disturbance equations are formulated as generally 
as possible within the framework of assumptions specified above. 
Elimination of other terms by examining the magnitudes of their 


dimensionless coefficients will lead to additional simplifications. 


2.3.1 Disturbance Equations 
With the assumptions discusscd in section 2.3 the disturbance 


equations for this problem become (dimensionally) 





Continuity: 


O22 
Homentun: 
TEC UEM Cic A SE 
" OU ne E UU OU M N 
ee e (2.23) 
A 5 2a AUS": SS 
217 da, y (0% =) 21$ 4 
t dy VAS ES) pd 2 E X 
o IO EE Wy OV eG cos XK - OP, 
JE uU ei az) Si o) 
ala, ww] , 31, 9% =/08,0v 95) (2.24) 
PAS + y] = al E 
+ dū An. dë An , 2 [32V oÑ 
du A AS Séi y EC E 
Y Ox Ay 0% va = A 
ENEE = "S 
= A du aw LU ow) pe eek 
ESE Oz 
PE Nc SU IE > |j dw + 28) (225) 
P5335: oy dy / 9? by 
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Energy: 
TE P ee, 
e E D SE T w 
=, TR, ad: WR). zZ (day, dw) 
Te al EISEN 
A E ON O dora A 
` FACES re 
(+ Te (RAÍZ) 
O O Ay (2.26) 


Notice that by virtue of the parallel flow assumptions, the 
mean flow coefficients of all disturbance terms are functions of 
y only, so these equations are now separable in terms of solutions 
of the form 
A id (x* Cos © FZ e c* 4^) 

Q) ge 
Noting that the assumed disturbance form is complex (with complex 
amplitude function a while the disturbance equations (2.22-2.26) 
are not, physical quantities can later be recovered by taking the 
real part of the complex solution. Thus, the disturbances are 
assumed to be sinusoidal waves of frequency w* = eee propagating 
at an angle 0 to the x axis. For this spatial analysis, E is 


* 
assumed real and A complex to admit solutions which are temporally 


k * 
only oscillatory but spatially amplified (et « 0), neutral (A, 0), 
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X 
or damped (e. >20): 


Such a formulation will reduce the set of disturbance 
equations from partial differential equations to ordinary ones in 


we 


y. Thus introducing 
f = Ue f(y) 
= Uri de) 
Y 
= We hy) 
e» Um Cy) 


£> 


x 


E 
T 


(2727 
^. ` 2 3 i exo | is (X ces 0 
E go U mG) | er 
ÓN LZ aen CG P 
Te (I sales) Gy) 
A = pe miy) 
e = Ce c On 
Ay 
ks Lia 
vhere 
bie Ue (x) coSO + We sın © (2.28) 


and defining for the mean and disturbance velocities parallel to 


the surface, respectively, 


+. Le) e kri Ier ef 
| + tanWtan® EE 


and 
u(y) + W (y) tan W tan O 


WC | « tan V tan O SE 





W 
where tan Y = 





and 9 is the angle between the x axis and the 
e 


direction of wave propagation, then the disturbance equations 





. become: 
Continuity 
ir (We) E LJ. i (50) = O 
C (2-915 
em nern EN 
SE + (Ed) Ww" = = Sex t sinX cos © 
E m p | Fz PES ( zuel 
+ , E M PA (Fr ec y D 
—— m => 
e [se x K e 
and E 
2 E V NM PE _ / 
La (FOX Wee) = eet em 
auc i NT E mw” 
Re gi gt E Re id 


(21323) 


A 
All 2} 
| 
Dee 
3 
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i 
Q 
+ 
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+ 
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Energy 


ges 


Cp | (weit + EBH" | - WEST en 


‘zee @| (lo geb cl W| + zen a (ta fl, ` (le) 


E 7 "ses [e en] - ge] 





where differentiation is again with respect to n to permit in- 
formation obtained from the mean equations to be used directly in 
the numerical solution of the disturbance equations, and vhere 


all fluid property variation is of the form 


A = = ex cab +Zsınd - c*£) 


CO 


"e 
~ 


The dimensionless wave number and Reynolds number are defined by 








E Sg 
ds 
and 
E ed 
cas v ; 
co Se 


Note that the composite x-z momentum equation is obtained by 
multiplying the transformed x and 2 momentum equations by cos@ and 


sind, and then adding. Furthermore, dilatory expansion energy 





drops out of the energy equation with linearization, and so must 
be considered a second order effect for tnis problem. 

To simplify these equations further, the assumptions made in 
section 2.2.1 regarding buoyancy, dissipation and expansion energy 
are used again, with some interesting implications. 

(1) Neglect of the dissipation energy terms (for =p. O) 
permits the energy equation to be written as a function of E and 
W rather than f, h, u, and V, thereby reducing the order of the 
system of equations by two from eight to six. In other words, when 
viscous dissipation is negligible in the disturbance energy equation 
then "the stability of a three-dimensional boundary layer to a 
plane-wave disturbance of arbitrary orientation reduces to a 
two-dimensional stability problem governed by the boundary-layer 
profile in the direction of wave propagation and by the mean tem- 
perature profile"[15]. Thus, alterations of the external flow due 
to changes in pressure gradient or cross flow would be manifest 
as chanpes to the variable coefficients in the disturbance equations. 

(2) When the expansion energy is negligible GERE ««1), 
variation of the thermodynamic pressure appears only in the momentum 
equations. Under these circumstances, explicit inclusion of the 
second viscosity coefficient, DE in the momentum equations can be 
avoided by using the "mechanical" rather than thermodynamic pressure. 
Such a change would not affect the mean pressure, for using 
equation (2.5) and the "parallel-flow'" and linearity assumptions, 


P, = P; however, it would result in a difference between the two 
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— m ~ 24 au, av , dw 
finetuation pressures KS PE ( DE UK SS = dy J a) 40 


or, in terms of equations (2,27) 


A ea 


As noted before, for the large Reynolds and small wave numbers 
anticipated, this pressure difference could only be significant 
for large values of X/u*. 

On the other hand, if the expansion energy is not negligible 
(say, (B T DE~OC1)), the appearance of X can still be avoided in 
the momentum equations. Elimination of the thermodynamic pressure 
fluctuations between the two momentum equations ( 2,33) 
yields a vorticity equation - an Orr-Sommerfeld type equation 
modified by the inclusion of fluid property terms - in which terms 
containing A identically cancel [16]. However, this procedure 
would require evaluation of higher order temperature derivatives 
of the fluid properties. Since the equations will ultimately be 
reduced to a system of first order equations for solution, such a 
step should be avoided if possible. 

In short, the second viscosity coefficient should influence 
only the disturbance thermodynamic pressure distribution. For 
water flow, it has already been noted that E must be very much less 
than unity, so the first procedure for elininating À is utilized 


for this analysis. 


(3) Buoyancy effects are negligible CR seco 
a re 


= 


UA 





Since the assumption was made in solving for the mean equations 
that buoyancy terms of order Gr/Re << 1, it is certainly reasonable 
to assume that the body force influence with terms now of order 
(Gr/& Re^) is even less significant. If buoyant effects were to be 
included, their contribution would appear indirectly through modifi- 
cations to the mean velocity and temperature profiles. Although 
it is known that for flovs driven primarily by body forces, heating 
from above stabilizes while heating from below destabilizes the 
flow throuch the influence of thermal instabilities, this effect 
will not be considered herein. 

Thus with the assumptions discussed above, the disturbance 


equations finally become: 


e . TEN. 
Continuity: ir(W-c) + ou ES He Eg) =O (2537) 
L 


Moren tum: 


| J Se 
E Brem, Eo] pl gm W by sop 





(2739) 





an een 
ENS ET Së WË Ee (2.40) 
o Ke Fr e C 
Since these equations are not in a form suitable for 
numerical integration, they are rearranged and recast as a system 
of six, first-order equations.  Denoting the dependent variables 


as 


IN 


= e js "T 
Eu 
T. 


(2.41) 


— 
o 


N N 


6 


E 

Bor 

É ds 
Si 

the six equations can be written as 


6 
/ 
a =) Cr Se (2.42) 
Ei 


where E is a 6x6 coefficient matrix with the following nonzero 


terms: 


T2 


ie 


mas wo 
C = = Y 
d € 
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Z', is obtained from the composite x-z momentum ecuation; 


Z from the continuity equation; 2 , from the energy equation; 


t 
ros 


and Zi 


, from the y momentum equation. 


2.3.2 Disturbance Boundary Conditions 

The disturbance velocity boundary conditions used in this 
analysis are easily obtained by requiring that there be no relative 
motion between the fluid and the solid surface et this interface 
(no slip), and that the velocity normal to the surface must also 


vanish there (impermeability). In notational form, at y = 0, 


w-Q (no slip) 


nu Y 
vz=0 (impermeabili ty) 


TANDO CA 


li 
e 


or at n 


The thermal boundary condition at the wall requires more 
careful consideration (see also discussion by Dunn and Lin [17]) 
The scaling procedure used in Appendix A indicates that for a 
region very close to the surface (y/ó << 1), the energy equation 
reduces to the one-dimensional, unsteady conduction equation 
(dimensionally) 

EE IE dT) 
Me a 2 1ES + 
O 
The inportant point to note is that longitudinal temperature 
fluctuations are negligible in this region. Thus, solving this 


equation in the solid and requiring that disturbance temperatures 
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and heat fluxes be continvous at the liquid - solid interface 


(y = 0), the thermal boundary condition becomes 





K — GE 4 rip- 
GE De e $i i 
W 





where 
oe and the subscript "s" 
2 ges EU denotes properties of the 
(^s e o? solid. 


It is simple enough analytically to require that the thermal 
inertial of the solid be sufficiently large that thermal fluctua- 
tions die out in the solid very close to the surface. Under these 
conditions, E BE plane boundary condition would be 1(0) = O 
(i.e., ez = eU» ). Realistically, however, the full 
boundary condition should be closely scrutinized when applying 
this analysis to a particular solid, at a specified frequency, etc. 


Yor this analysis, the condition 


To) = 0 
(2.44) 
is chosen. 


Further it is assumed that all disturbances are bounded at 


large distances from the surface = 
^M ne me y - X 
geo U, I D 


ME (2.45) 


) 





2.4 Disturbance Field Energy Balance 

Although mathematically the stability characteristics of 
the disturbance flow field can be completely described (within 
the framework of specified assumptions) by the equations formulated 
in the three preceding sections, it is instructive also te develop 
the theory on a more physical foundation as well. Insight into the 
physical mechanisms at work in the boundary layer is thus sought by 
exanining the production, dissipation and transfer of the dis- 
turbance energies. 

The appropriate equations are obtained by first multiplying 
each of the x, y, and z disturbance momentum equations (2.23-2,25) 
by its corresponding disturbance velocity, u, v, and w, and then 
summing the three. After some algebraic manipulation, a meaningful 


equation is obtained and expressed (dimensionally) as: 











(1) aV 
nn? ry 2 yz y ü ^J mg W 
O /U 4*V 4 W ie 6 (UV tM 
DLE ( ==: y dy 
(3) (4) 
SS ME / 
3 fe uU = o W > Wes de 
— (ën TL SR - BD +4 ` 
Z Y X 
(5) e 
Eure A Te 
Mc ~ oJ = ru "T ~ a T Noc ru 
+ Qe AStA- A) Eo. 
dy l OX oz X B I f (2.46) 
(7) (8) 
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where v=o 
M NI E 
D T DE W 53 
= E + E è = (ee 2 x 
z M a X ay 3 
es ow +20) t. S Ou GE 
2 oz ? Oz OX 
= (ev - 98) 
E òX O y 


Assumine only tuo-dimensional mean and disturbance flows 
c 3 


equation (2.46) reduces to: 
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The time-averaged version of equation (2.47) written nondimension- 


ally in terms of disturbance amplitude functions (see equations 
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(2.48) 











and * denotes the complex conjugate. 


It is instructive to examine each term at this point to 
identify the physical process it represents and assess its relative 
importance. Thus, the indicated terms represent: 

(1) the rate of change in kinetic energy of the disturbance 


field. Note that for a neutral disturbance, the averaged value 


is zero so there is no net energy chanfe in the boundary layer. 





(2) the energy drawn from the mean motion by the working of 
Reynolds stresses against the gradients of the mean velocity. As 
this is the only source of energy for the disturbance field of 
the unheated flow, it can also be expected to be significant when 
there is heat transfer. Furthermore, it is expected to be greatest 
and positive in the region roughly between the wall and the 
critical point, signifying energy transfer (Sion ne mean to the 
disturbance flow, but small end cuite possibly negative from there 
to the boundary layer edge, 

(3) the energy extracted from the mean flow due to the 
working of the disturbance viscous stresses against the mean 
velocity gradients. by and large, it is expected to have the same 
sign as term (2) and so should represent. primarily a source of 
energy for the disturbance field. It should be most significant 
in regions where temperature fluctuations are largest since the 
viscosity and temperature fluctuations are proportional. 

(4) the rate at which the pressure fluctuations do work. 
Physically this energy is produced by two different mechanisms: 

(a) the periodic expansion and dilation of each fluid element 
(caused by density fluctuations) in the presence of a disturbance 
pressure field, and (b) the movement of this same fluid element 
along a mean density gradient also in the presence of the dis- 
turbance pressure. Note from equation (2.48) that for a neutral 
disturbance, the contribution from the first mechanism will vanish 


at both the wall and the critical] layer. 





ao 


(5) change in energy resulting from a nonzero mean viscosity 
gradient. As may be seen from equation (2.48), this term is of 


4 ot, 


order times that of the major production terms (2). 


Nel 
Since the general order of e decreases and Re increases with 
moderate heating, it is expected that this term should be of 
negligible significance. 

(6) mechanical energy lost from the disturbance field by 
viscous friction. This loss is manifest physically as a heat 
source which in turn produces minute changes in the disturbance 
temperature (and so density) level. Since the term is negative 
definite, it will alvays represent a dissipation and should be most 
important near the vall. 

(7) energy dissipation from thermal fluctuations.  Specifi- 
cally, this is the energy lost associated vith the dilation of the 
fluid elements due to their expansion and contraction with changes 
in density (caused by changes in temperature). This effect should 
be small in liquid boundary layers where the amplitude of dis- 
turbance density fluctuations is relatively small, even with 
significant mean thermal gradients. 

(8) rate of transfer of pressure and dilatory energy, 
vorticity, and disturbance viscous shear stress. If the disturbance 
field is neutral (x, = 0), no energy is transferred parallel to the 
plate. Further, if this term is integrated from the vaJl to some 
large distance (n > 9), the remaining energy contributions will 


vanish identically due to the homogeneous boundary conditions. Thus, 
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the term acts only to transfer energy from one position in the 
boundary layer to another, and for a neutral disturbance makes no 
contribution to the net energy. 

Two examples are now examined to more clearly identify the 
terms of interest in this investigation. First, equation (2.48) 
is reduced to a form compatible with the assumptions made by 
Wazzan, et.al; namely in the disturbance equations, p = constant 

e = ft op 

(but variable in the mean flow), u = uy (y), and 1 = 0 (so p=) = OQ). 
Thus, for a neutral disturbance, the y-integrated, nondimensional, 


energy equation is composed of only the production and dissipation 


contributions: 





(2) Ante m 
O = - 32 Re a (204) > 4p (gd, ), 
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(2.49) 


(Note that the variable p still appears in this equation since 


d = z 
Ær erentiation is with respect to n; i.e. daer) Ts’ Q Sch 


Secondly, assuning for this analysis all mean and disturbance 





quantities may vary, the equivalent equation becomes 
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where “X (q) is here slightly simplified to 
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CHAPTER III 


NUMERICAL SOLUTION OF THE STABILITY PROBLEM 


J.1 Thermodynamic Properties and Transport Coerficients 

Before numerical boundary layer integration can be performed 
for either the mean or disturbance equations, nc for the 
thermodynamic and transport properties of the liquid and their 
required derivatives must be specified at eacb step. Consistent 
with assumptions made thus far that the property variations are 
independent of pressure (also for the boundary layer approximation, 
the mean pressure is constant in the boundary layer anyway), these 
quantities can depend on position only through values of the mean 
temperature. For the mean flow, such temperature dependence 
couples the momentum and energy equations. More specifically, for 
the stability problem, accurate information is needed not only for 
the variable properties of the mean flow, but for the disturbance 
quantities r, m, and x as well (see equations (2.20 and 2.35)). 

Since the solution process for the mean flow equations 
requires changes of the temperature at a particular physical 
position both during the iterative search for unknovn wall boundary 
conditions (discussed in section 3.2) and later with the choice of 
step size for the converged solution, discrete input to a numerical 
scheme of required fluid properties and derivatives is impossible. 


That is to say, this information is available in general ouly by 
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interpolating tabulated experimental data at each step. It 


ollows, then, that the temperature derivatives must be obtained by 


En 


numerical differentiation of the discrete set of data points, and 
so must be interpclated as vell. 

Water is chosen as the representative liquid for numerical 
purposes to maximize applicability of the results of this investi- 
gation (in vhich thernal fluctuations are considered) to situations 
occurring in nature, but primarily to enable comparison with the 
results of the fourth order stability analysis of Wazzan, et.al, 
[4]. In general to obtain the "best" information available, 
values for ET D ,u , and k are taken from recently published 
measurements and/or evaluations of preexisting data. Temperature 
derivatives of these propertics are obtained directly by 
analytically differentiating the empirical property~temperature 
relationship proposed for each. Note that this procedure is in 
contradistinction to first estimating derivatives directly from 
the data (with the method of cubic splines, for example) and then 
smoothing any spurious fluctuations. 

A more extensive discussion of the sources and accuracy of 
input property data, alternate means considered for curve-fitting 
tabulated data and obtaining high quality derivatives, and graphical 
presentation of the properties and temperature derivatives for 


water are given in Appendix B. 





n 


3.2 Mean Equations 
The equations describing the mean two-dimensional, boundary- 
layer flow over a heated, flat plate are obtained directly from 


equations (2220 eas 


H 
(52 E") "ze EF"=0 
and. (3.1) 
i TE A / 
(ek H $ Gt pedo 
A 
with boundary conditions specified by equations (2.21) 
/ 
"e E | 
and (2:2) 


The numerical solution to these equations introduced several 
difficulties characteristic of a boundary value problem, in 
addition to those associated with SES nonlinearity and coupling 
through fluid property variation. The method of Nachtsheim and 
Swigert [18] adopted for this analysis is found to satisfactorily 
deal wich ai difficulties, especially: determination of where 

to stop the integration for the unconverged solution, and 
approximation and refinement of initial conditions. 


First, to obtain a unique solution of these equations, 





external conditions must be "satisfied" for an undetermined, but 


y 


finite n =n > The Asymptotic convergence (where, F' > 1 and 


max 


n ) is assured 


H>+0as n > becomes F'  andH 0 at n -— 


Pep ly by requiring velocity and temperature gradients vanish at 
the "truncated infinity" (F" = H' = 0 at n= pc Note from 
equations (3.1) that these additional conditions cause all higher 
derivatives to vanish there as well. However, the value of [E 
(or equivalently, that n beyond which no appreciable change is 
observed in the solution) still remains unknown. If the a priori 
specification of Ee is too small, boundary laver integration is 
terminated toe soon and boundary conditions cannot be satisfied; 
at RE is chosen too large, computation time becomes excessive. 
Secondly, numerical integration of the mean flow equations 
requires specification of as many conditions at the wall as there 
are conditions to be satisfied at the boundary layer edge. Since 
the required additional wall conditions (F and H^) are initially 


unknown, they must be first estimated and then iteratively adjusted 


until the asymptotic conditions can be satisfied, or 


EAM f £ 
a (Fy , im : n) : 


ait bi = 3, 
Bue B los E (3.3) 
with the supplemental conditions 


ESOS = 0 
Fi W eeneg 


1 fp iR E 
: Po is D) H 


In other words, the system of equations is solved as an initial 
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value problem with the wall gradients as the specified initial 
values. Due to the extreme sensitivity of equations (3.1) to 
values of De and Ho even falr estimates of these parameters are 
Biskely to cause the numerical boundary layer integration to "blow 
up" prior to reaching the prescribed value of fl RE Thus no means 
would be available for refining initial estimates of the wall 
gradients. 

The mean equations are thus solved in the following fashion. 
With an estimate of initial conditions, 2 and n^ Runge-Kutta and 
Adams-Moulton integration techniques are applied over an abbreviated 
n range. At the end of this range, the initial conditions are 
refined by determining the least-square solution for discrepancies 
between the computed variables and the proper asymptotic values. 
Wnen successive iterations over this interval fail to appreciably 
change the initial conditions, an error estimate E, the sum of the 
squares of the deviation of computed quantities from their 
asymptotic values, is evaluated: 
2 


2 2 


Baal = pty? +H” + FU + Hi (3.4) 


Hf this error is sufficiently small, then the solution is considered 
to be found. If such is not the case, then the n range is increased 
and the above procedure repeated. The boundary layer is thus 
traversed in a stepwise fashion until errors are within acceptable 
limits. The last n range then determines Be automatically. 


In summary, this method yielded relatively rapid, unique 


convergence to the proper solution and seemed to be insensitive to 
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initial estimates of the wall gradients. 

Once solution of the mean flow equations is accomplished, 
the variable coefficients for the disturbance equations can be 
E puted E ect. and the latter solved in the manner described in 
the next section. However, other incidental calculations can 
easily be performed as well. For example, to find ns for the 
unheated case (or more specifically, when density is considered 


constant), from equation (2.18) 


co Imax , 
wl (Uc uA = ae ao man = F Oman) 


For the heated case, though, the solution is not immediately 


availeble from variables computed in the boundary Jayer integration: 


e SS J SE? u FO nox) (3.4) 


Instead, using stored values of p (n) and the composite Simpson's 
Rule formula [19, p.79], the required € em is determined 
numerically.  Computed values of n^ are tabulated in Table E.3 and 
plotted in Figure C.2. An analysis in section E.] indicates an 


«bs - 10 to be adequate for the boundary layer integration of the 


mean equations. 





Conversion back to physical coordinates must be similarly 


computed. From equations (2.16a) and (2.17) 
5 | 
kc | z (3.5) 
* . 
(CH 


(3.6) 


so that 5 _ El CR 
5” de SS 


a constant, once the mean flow equations are solved and the 
* 


boundary layer edge is specified. A plot of the kin and 3. 
Ò 


relationships for the mean flows analyzed is given in Figure Cl, 


3.3 Disturbance Equations | 

The equations describing the hydrodynamic spatial stability 
of two-dimensional disturbances in a two-dimensional, heated, 
liquid boundary layer are obtained from equations (2.37-2.40) by 
setting 6 = 0 and y = O (see equations (2.29 and 2.30)), and 


requiring that Re and w = «c be real and & be complex. Thus, 





rol 


Continuity: 
[o = / 
Momentum: 
if (G-c) +@a)a = E E = por +3] — 
zur 
J 4 ee E ; O 
Ø Ee E 
| =a {ene pif: 2 f (3.6) 
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Energy: 


= 
| i 1 NE m / N, 
e [i-e * eóH] siege eet] SS j 


which satisfy a set. of homogenous boundary conditions obtained from 


puusrions (2.052.605) 


and (3.7) 


B Jonn 


Since these equations and boundary conditions are insufficient to 
establish a nontrivial solution, the system is cast os an 
eigenvalue problem (i.e., nonzero solutions satisfying boundary 
conditions are obtainable only for selected combinations of the 
parameters a, w, and Re). 

Since the complexity of the equations precludes exact 
analytical solutions within the boundary layer, numerical techniques 
are applied. Of numerous possible methods previously utilized for 
solving this type of problem, the orthonormalized, differential 
method used for this investigation seems to be one of the most 
efficient and general [20]. Schematically, the employed scheme has 


the following steps: 





(1) estimate an eigenvalue; 
(2) starting from solutions valid in external region, 
integrate numerically to the wall and try to satisfy 


the specified boundary conditions there as well; 


(3) refine the initial eigenvalue estimate and reinte- 
erate as Jn (2); 


The final two steps are repeated until "satisfaction" of wall 
boundary conditions is attained. Although this procedure is 
conceptually quite straightforward, its implementation is not. The 
ensuing discussion elaborates on the difficulties associated with 
each of the steps enumerated above. 

As noted previously, there sre four parameters associated 


with this eigenvalue problem, X , «&,, o , and Re, In subsequent 
R 


i 
discussion, "the eigenvalue" shouid be interpreted as those two 
real parameters which are iteratively adjusted to satisfy the real 
and imaginary parts of a wall boundary condition. The choice of 
the two fixed parameters is determined by information sought from 
the analysis; for example, to obtain curves of constant amplifi- 
cation, it is found convenient to fix Ay and Re,and then search 
for the proper values of a and o0 . Although inclusion of a 
fluctuating temperature field for this investigation is expected 
to somevhat modify the existing results of Mazzan, Okamura, and 
Smith [4] (in which no such field was considered), the initial 
eigenvalue estimate for a specified T and T. is taken from the 


latter. Once several "exact" eigenvalues are determined, these 


can be extrapolated to estimate others in the same neighborhood. 





problem is row converted to a boundary value problem. 


Sl 


With all pertinent parameters specified, the eigenvalue 


Pause 


remains to mathematically describe the asymptotic condition 


specified for the boundary layer edge. 


Since the mean flow co- 


efficients in the disturbance equations are all assumed to be 


constant outside the boundary layer, an exact solution can be 


Specifically, equations (2.42) reduce to 


obtained for n » Dat: 


Z - -i é, - 16:3 e 


2 
Zz =z, 
; E = 
Z= -|ia (i-o) P EE 
AU : 
5"& (8 7 94 
Zs = [ iaRe(j-c) + 97] Z, - a" (48) 


A + ¡ale P- (1- c)|£, 


6 


E 
or after some manipulation 


uil m t jl aw? es 
2 ES EE cu. + ee) 


and 


2 = e = O 


(3. 8) 


(G-Z -wz | = © 


(3.9) 


(3.10) 





where 
DM NM ES 
2 Se Ga, (3.11) 
O zu + ¿el Re te (4-e) 


From Eqn. (3.10) it is easily shown that to satisfy the external 


boundary conditions (equation (3.7) ), for the JER solution vector 


Ec Ober en 


for n» n, ond real B » U 
(2512) 


Thus equation (3.9) is an ordinary, nonhomogeneous differential 
equation with constant coefficients which when satisfying the same 
conditions of boundedness has the solution 
e d 
cO (GO — Ga: SE io (SB > 
ZU. C. E Sa e 


2 


7 
(Ea) 


(3:13) 
for n > no and real a, y > 0. 
From equation (3.8) and (3.13) 
Z0 ig Ce e MON 
2 (9.14) 


et 


mor n > ne and real a, y > 0. 
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The values for Lys TE and Ze can be expresscd in terms of Zp bos 





pos O LE (ia TG 


Sëch m: z4 - (a4 y)Z9 [eerie Gee 


So) 


As indicated, the three acceptable, linearly independent 
fundamental solutions E , e E Ene , and coe 
can be combined to obtain a generali solution to the simplificd 
disturbance equations outside the boundary layer. 

With an eigenvalue estimate and specified values for the 
y) 


three linearly independent solution vectors at the boundary 


q = 6 


layer edge (obtained from eqns. (3.12-3.15) by setting C ij 


1, j = 1,2,3), a fourth order Runge-Kutta method is implemented 

to numerically integrate simultaneously the disturbance cquation 
solution vectors to the wall, n = 0. After each integration step, 
however, the original linear independence is found to be so greatly 
diminished that the solution vectors become computationally 
dependent before reaching the wall. ‘Thus, insufficient information 


is available there to satisfy boundary conditions. 





In complex vector space, this can be viewed as a rotation of one or 
more cf the solution vectors from an orthogonal to a colinear 
configuration, or a reduction of the included angle from ninety to 
zero degrees. The difficulty is explained as follows. Although 
the exponentially growing solutions to eguations (3.16) can be 
excluded from the boundary conditions in a mathematically exact 
Way, their influence cannot be eliminated from the numerical 
integration inside the boundary layer. Here, truncation or round- 
off errors provide suitable initial conditions for subsequent 
integration steps. As the computer will allow tne most general 
solution, the rapidly growing portion eventually dominates the 
three solution vectors (parasitic error). To avoid this difficulty, 
an orthoromal basis is reestablished whenever complete loss of 
linear independence is imminent. Gersting and Jankowski [20,21] 
contrast this '"near-orthonormalized integration" to the "fully- 
orthonormalized integration" used by Vazzan, et. al. [3], in which 
solution vectors are kept orthonormal at each step. 

Integration thus proceeds in the folloving manner. The 
three basic solution vectors z) are integrated using a Runge- 
Kutta scheme with a fixed step size from the boundary layer edge 
tovard the wallimtil, at mesh point n4» the attendent par3sitic 


error reduces the angle between any two solution vectors below a 


Epceificd lizit, or 


» E E DROP ios 
iD m) al de (3.20) 


min COs 





where Gil, z^ revresents the complex inner product of the 

. vectors i and j. There, the Gram-Schmidt algorithm for 
orthonormalizing a set of vectors ís applied [22]. When 
eigenfunctions are to be recovered, tne orthonormalizing matrix, 
available as a nonessential by-product of the orthonormalization 
process, is also computed and stored for later use. The integra- 
tior is then continued from ny using as initial conditions the 
orthonormalized basis vectors, zu until at mesh point Nos the 
angle between any two solution vectors is again less than 2. This 
procedure is repeated until the vall is reached, at which point 

a final orthonormalization is performed. 


An evaluation of the proper step size, h,, and orthonormali- 


d^ 
zation angle criterion, &, required for this integration is given 
in sections E.2 and E.3 of Appendix E respectively. Based on that 


evaluation, the values of h .02 and 2 = 45° were used to generate 


nes 
the results presented in Chapter IV. (Since the integration scheme 
used requires computation of derivatives at half-step intervals, 

the mean equations thus must be integrated with step size DM - 

1/2 hy. 

At the wall, the three independent solution vectors are 
linearly combiued in an attempt to satisfy the homogeneous boundary 
conditions, equations (3.7), for that particular eigenvalue estimate. 
In general, such satisfaction is not achieved immediately unless the 


eigenvalue is "exact," and so some iterative correction scheme must 


be implemented. Thus, for the vector components 





EO. e) 
le E; o EV) 


it is required that (see also [23]) 


Z, (0) + B, Z (o =O 


and (3218) 


XE NO 


= 
where specification of the arbitrary constant B, = u (0) is found 


1 
to give good eigenvalue convergence using the search procedures 
discussed below. Following Wazzan, et. al. [23], the quantity 
actually used to clearly distinguish an eigensolution from all 
others is the complex, surface-normal, boundary layer admittance 


t (w, a, Re) = 


(o) (3.19) 


for, the smaller the "zero" components of the eigensolution 

(2, (0), 2,0), and 2,(0)) are than the nonzero ones (2, (0),2. (0), and 
Ze fei), the more accurate will be the eigenvalue [24, p. 280]. 
Requiring both the real and imaginary parts of Y. vanish establishes 


tvo real equations in- the two adjustable parameters (or the complex 





"eigenvalue"). 

For iterative refinement of the eigenvalues (i.e., finding 
the zero(s) of Y)» two schemes were employed: the plane fit of 
‘Wazzan, et.al [23], and the quadratic complex curve fit of 


Muller [25]. In the former, a plane of the form 


"eigenvalue" = a, + a Y, ^ aX. 


l oR 2 oi 

is fitted to three successive eigenvalue estimates and their 
corresponding computed aamittances. The corrected eigenvalue 
estimate for the next iteration is then simply E When H is 
sufficiently small, Muller's method is used to expedite eigen- 
value convergence. This scheme fits a quadratic curve in complex 


space, also through the last three eigenvalue estimates and com- 


puted admittances, of the form 


Y oz b + b, (eigenvalue) + b, Cena 


0 I 


Setting D = 0, this equation can be inverted using the standard 
quadratic formula to give an eigenvalue estimate for the next 
iteration. Note that both techniques have the advantage that no 
derívative of Lë and only one evaluation of P is required per 
iteration (in contrast, say, to Newton's method). A complete 
discussion of the more subtle applícations of these tuo procedures 
r an Efficient eigenvalue search is relegated to section E.4. 
When the “eigenvalue” is determined to within desired 
accuracy, the eigenfunctions can be recovered using the stored 
Srthonormalizing matrix and solution vectors [22]. For standardi- 


zation, the magnitude of these eigenfunctions are adjusted so as to 





make the amplitude of the pressure fluctuation equal to unity at 
the vall. If desired, energy production, dissipation and transfer 


terms can then be evaluated directly (see section 2.4). 





CHAPTER IV 


RESULTS 


The hydrodynamic stability of the boundary layer formed on 
a heated, constant-temperature flat plate (M = 0) in two-dimensional 
water flow and with two-dimensional disturbances (Uv = 0 and 6 = 0) 
is examined and the results presented in this section. Specific 
comparison is souzht with the results of Wazzan, et.al [4] who 
analyzed this same configuration, but who assumed only isothermal 
disturbances of the mean variables. 

Rather than regenerating in this investigation the Stability 
characteristics for the entire range of wall temperatures considered 
by Wazzan, et.al., selective checks are made to compare the results 
of the two. Accordingly, since those authors found that increased 
surface heating at first stabilized and then destabilized the 
boundary layer flow (i.e., peaking of Re - that value of 


Mia Crit, 


Re below which disturbances are damped for all wave numbers and 
frequencies), with the extremums occurring at about T. = 140 F 

T. -= 60°F), the numericsl calculations herein will be limited 
primarily to wall temperatures within the normal liquid range of 
water and including that temperature ranye where the peaking occurs: 
specifically, T,= OG toe wand AO Tr nih T = 60°F When 


detailed characteristic behavior is to be illustrated, the case of 


L, = 90°F will be used, as this is conceivably within the ranpe of 
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temperatures for which experimentai comparison can be made, 

Àn accurate solution to the mean fiov equations is essential 
for solving the disturbance equations. The procedure and results 
‘for the mean flow solution are described in Section 3.2 and 
Appendix C. In the latter both detailed boundary layer distri- 
butions of all the mean properties and their derivatives required 
for the coefficients of the disturbance equations, and variations 
of surface velocity and temperature gradients witn the wall 
temperature are graphically presented. 

The mean velocity and temperature boundary layer profiles to 
be examined in this section are repeated in Figure 4.1. Note that 
the temperature boundary layer is much thinner than that of the 
velocity boundary layer so thermal influences are restricted to 
repions close to the wall. This observation is related to the later 
understanding of the relative importance of thermal disturbances in 
boundary layer stability calculations. 

Keeping in mind the role that a linear stability analysis 
should play in only suggesting those conditions which might enhance 
or delay the onset of laminar-turbulent transition (there is no 
definitive link between the point of instability predicted by such 
En amalvsis and the point of transition), this investigation is 
restricted to regions which could conceivabiy still be laminar. 
Thus, although determination of stability characteristics for 
excessively large Re might be an interesting numerical or mathemati- 


cal cxercise, the results are not of interest if transition has 
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occurred. 

Much of tnis investigation was conducted neglecting the 
effects of a fluctuating density field (r = 0); that is, only 
fluctuations of viscosity m(y) and thermal conductivity k(y) were 
included in the disturbance momentum and energy equations, respec- 
tively. Such an assumption seemed reasonable due to the relative 
constancy of the density of water over the normal liquid range. 
Subsequent calculations with variable density nevertheless were 
made and their results are included. 

Unless otherwise indicated, the property-temperature varia- 
tion for water is that specified in Appendix B.2. Although this 
variation differed slightly from that of Kau»s and Smith [7] (used 
by Wazzan, et.al. in their analysis) as shown in Figure D.7, there 
is a marked effect on some of the calculated stability character- 


istics, as wiil be shown presently. 


A A mm e ` mm 


As suggested in previous discussion, for a specified set 
of mean flow parameters (T I and M), the eigensolutions to the 
homogeneous disturbance equations and their boundary conditions 
describe a surface in the four-dimensional parameter space 
(e, >» Gp» Us $ and Rc). It is of interest to locate minina and 
maxima on this surface and to trace the path followed by an actual 
disturbance as it propapates through the boundary layer (e.g., to 


ascertain whether an amplification repion (a, < 0) is traversed, 





and if so, the rate of amplification expected). The stability 
characteristics associated with neutral and amplified disturbances 


are discussed below. 


4.1.1 Neutral Disturbances (ay = 0) 


e mm e —À re E a a e Re m 


Perhaps the most obvious quantitative puide to the effect of 
heating on boundary layer stabilization is the amount vith which 
neutral curves are shiíted with changes in the wall-to-free stream 
temperature difference. The movement of these stability loops is 
characterized by changes in the minimum critical Reynolds number, 


Re . „the maximum wave number a , and the maximum fre- 
Dan.crit, RM 


quency, w_ (above which all disturbances are damped). As may 
max 

be seen from Figures 4.2 and 4.3, with only moderate heating, both 

the fourth order system of equations of Wazzan, et.al. [4] and the 


sixth order system of this investigation (without density fluctua- 


tions, however) predict a significant delay in the onset of 


) 


b da lay instabilit vide d increased Re . 
oundary ayer in ility (evidenced by the increasec AA 


and a reduced range of frequencies and wave numbers which could 

become unstable. As the wall-to-free stream temperature differcnce 

is further increased, however, both the fourth and sixth order 

systems consistently indicate a peaking and subsequent decrease of 

Re .. . Based on these calculations for water at the indicated 
EDU crPIt. 

temperature levels, the following conclusions may be drawn. 


(1) It would be impossible to achieve complete boundary 


layer stabilization by indefinitely increasing; the heat rate; 
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EL) While the calculations which lnelide thermal disturbances 


do predict slightly larger values of Re the quel cative 


EE 
variation of this auantity with P (Gor fixed a 60°F) is not 
altered. 

The values of Re are easily obtained for specified 


PAN. Crit.: 


mean flow parameters by fixing a, = 0 and a,, then adjusting esti- 
A 


R 
mates of Re and w, until boundary conditions en be satisfied. An 
example of this procedure is shown in Figure 4.4 and 4.5 for several 
different cases. 

For no other calculation is the influence of the chosen pro- 
perty-temperature relationship or assumed property fluctuations more 
pronounced than in this one. Tor example, from these two figures 


it may be seen that: 


(1) use of the property variation of Kauns and Smith [7] 
ei C S a shifil of he first above c = 90°F) 


Wi Cri, 
and then below (T 


H 
> 


200°F) that computed with the 
variation specified in Appendix B.2; 


(2) for either property variation, calculations including 
thermal disturbances (but with r = 0) predict a 


e , . roughly l10Z grcater than that computed 
Mun Crit pn E p 


using the fourth order eguation of Wazzan, et.el. 
(equation (1.2)); 


(3) inclusion of density fluctuations in the calculations 
os well, offsets the gains iu stability predicted above 
Ina as T is increased. 


An explanation for the relatively nronounced influence of 
density fluctuations on Re . calculations is found in the 
Miser ic. 


distribution of disturbance properties in the stability equations. 


Although the density fluctuations may be comparatively small (by 
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as much as several orders of magnitude less than the viscosity 
fluctuations - see Figures 4.9 and 4.10, for example), they do 
appear in both the inviscid and viscous parts of the stability 
equations, whereas viscosity Fluctuations are restricted only to 
the latter. As temperatures are increased, the disparity between 
the two disturbance quentities diminishes; eventually, disturbance 
density contributions dominate those from all other property 
Peictwations in the calculations. 

However, as shown in Figures 4.4 and 4.5, while the fluctu- 
ating viscosity terms tend to enhance stability, the density terms 
degrade it. It is important to note that over the normal liquid 
range of water, results from the sixth order system of equations 
(even including density disturbances) predict a slightly more stable 


boundary layer than does the fourth order system for which r (y) z 0. 


4.1.2 Amplified Disturbances (a, < 0) 


- 


The tendency to transition is also a function of amplifica- 
tion rates and so attention will now be piven to the amplified 
portion of the stability loop (a, < 0). 

A typical set of stability curves illustrating the topologi- 
cal nature of the eigenvalue variation for T = 9Q°r, ne: 60°F, and 
M= 0 (r = 0) is shown in Figure 4.6. These curves are similar in 
gencral shape to those for the unheated Blasius boundary layer 
although dramatically different ín magnitude (see reference [26, p.90] 


for example). 
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A much more significant indication of boundary layer stabi- 


lization by heating than Re is given by the maximum 


minerit: 
amplification rate, II established for specified mean flow 
i » 
conditions. The parametric location at which this maximum occurs 
is determined by a more involved procedure than that for finding 
e . meee as 1 biustroved in Ficure 4.7. More specifically, the 
nin. erit. 
amplification rates for specified Re are computed and individual 
local maxima determined.  Folloving the locus of these individual 
maxima, the largest value attained is designated the maximum ampli- 
fication rate and its associated wave number (and frequency) are 
obtained graphically. While the flatness of this locus curve in the 


neighborhood of -(a, /Re) allows relatively accurate evaluation 


max 
of that quantity, it does, however, preclude similar accuracy for 
De (and o, ) there. Physically this implies that there is effectively 
a band of wave numbers (and frequencies) which will receive nearly 
maximum amplification. To complete the eigenvalue set, Re can be 
determined from Figure 4.6 by again following the locus of indivi- 
dual maxima until the required value of Op is reached. 

Due to the numerical effort required to compute the maximun 
amplification rate, only one set of calculations was made for 
m - 90°F and D 200°F (T, = 60%P), using the presumably more 
accurate property-terperature relationships of Appendix B.2(assuming 
r: 0). Values obtained for -(a,/Re) ax x 10° of 26.0 [or the Former 


and 19.5 for the latter can be compared with the corresponding 


values tabulated by Wazzan, et.al. [4] of 28 and 20, and the valuc 
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for the unheated (Blasius) boundary of 740. 

Thus, not only is the onset of disturbance instability delayed 
and the range of disturbance wave numbers and frequencies receiving 
EE agoen reduced as pointed out in the preceeding section, but 
the amplification rates to which disturbances will be subjected are 
also dramatically reduced. These featurcs all point to a more 


stable boundary layer flow. 


4.2 Disturbance Amplitude and Energy Distributions 

Finding no major modifications to the stability characteris- 
tics of a heated, water boundary layer by including in the calcula- 
tions the effects of a disturbance temperature field, an explanation 
is sought from the boundary layer distribution of disturbance 
quantities and the energies produced or dissipated by them. 

Due tc the homogeneity of the equations and boundary con- 
ditions, solutions to the disturbance equations are independent of 
an absolute emplitude (ss section 2.3.1). Thus, as a comparison 
is sought between disturbance levels for both the same and different 
mean thermal boundary conditions, it is necessary to resort to some 
physical reasoning to determine a normalizing quantity that would 
remain virtually constant in the comparison. Since there is no 
streamwise pressure gradient (i.e., M - O0 for all Re), the mean 
pressure vithin the boundary layer is assured constant, and the 
disturbance pressure is nonzero and slowly varying throughout 


(see Firures 4.8 through 4.10, for example), the amplitude of the 
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Figure 4,8 Neutral disturbance amplitude profilas for the 
unheated (Blasius) boundary layer. 
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disturbance pressure at the wall, |w (G)|, was chosen to normalize 
all disturbance quantities. Thus all comparisons made between 
different eigensolutions in the remainder of this section are pre- 
dicated on this assumption. 

To further investigate the importance..of density fluctua- 
tions, they are included in all calculations made in this section. 

Typical amplitudes of disturbance quantities found in this 
analysis are presented in Figures 4.8 through 4.10 for wall tempera- 
tures T = 60%, 90°, and 200"F (T. - 60'F), respectively. For 
standardization, the eigenvalues associated with each set of eigen- 
functions are all located on the upper branch of their respective 
neutral stability curves at frequencics near that receiving maximum 
amplification. 

After studying these figures, the following observations are 
made. 

(1) As anticipated from the relatively large Pr for water 
(i.e., such that the thermal boundary layer thickness is much less 
than the velocity onc), in the temperature range considered, thermal 
effects are restricted to a region very close to the wall. Since 
thermal fluctuations must vanish at the wall, even within this 
region, they attain fairly large enplitudes only within a far narrow- 
er range. 

(2) By far the most dominant property disturbance amplitude 
is that of viscosity. However, for relatively small surface heating 


M = 90°F), even the thermal conductivity disturbance amplitudes is 


0 S| oS a 
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commensurate vith that of the normal velocity in some regions. 

For large heating v = 200°F), when density fluctuations are 
significant, they are seen to be larger in amplitude than both the 
thermal conductivity and normal velocity disturbances. 

(le al 


(3) With increased heating, the ratio of | £ | 


mas max 


increases significantly. Further, assuming that the disturbance 
wall pressure amplitude is relatively insensitive to changes in 


mewithin the range 60° = T = 200°F, | £ | 
W 


" increases monotonically 
N 


max 
with T while Ie el... exhibits a maximum. 

(4) Just as the stability characteristics are only slightly 
modified by the inclusion of thermal disturbances in the analysis, 
so too are the disturbance velocity and pressure amplitude pro- 
files. Even the disparity shown (except in the region where thermal 
disturbances obviously dominate) can be attributed in part to the 
difference between the neutral eigenvalues required to satisfy the 
fourth and sixth order systems. 

Once the eigenvalues and eigenfunctions are obtained by the 
procedure outlines in section 2.3.1, it is a simple matter to com 
pute the energy balance developed in section 2.4. However, for ease 
of analysis, only those terms contributing a net increase or de- 
crease of energy of a neutral disturbance are evaluated (i.e., all 
terms of equation (2.50)). The results of such calculations are 
presented in Figures 4.]] through 4.16. For each of the cases 
E = 60°, 90%, 200°F (T_ = 60°F), the eigenvalues selected from the 


upper and lower branches of the neutral stabílity curves (see 
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Figure 4.16 Energy distributions of a neutral disturbance in the 


heated boundary layer: 


t= 200°F and T = 60%F (upper branch). 
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Figure 4.1) again are for frequencies for which disturbances are 
subjected to approximately maximum amplification when traversing 
the loop. 

From these six figures, the following observations are made. 
Ihe "terms" referred to are those designated in equation (2.50). 


(1) The dominant terms in the neutral disturbance energy 
balance are the traditional ones: The Reynolds 
stress production (term (2)) in the "eus between 
the wall and the neighborhood of u = C,, balanced by 
the dissipation (term (6)) which is IE: Jarge 
near the vall. 


(2) Dissipation due to a nonuniform mean disturbance 
density (term (7)) is always negligible. Although 
not nearly as small, energy production due to a mean 
viscosity gradient (term (5)) is also insignificant. 
Neither are large enough to appear in Figures 4.11 - 
416, 


(3) The disturbance shear stress (term (3)), although 
dissipative in some regions, represents primarily 
energy production for the disturbance field. In 
all the cases examined, it produced roughly 10% of 
the energy, with the Reynolds stzess production 
(term (3)) supplying the rest. Recall that this 


is approximately the same amount that Re . 
IP o 


was increased for sixth order system calculations 
(r = 0) beyond that calculated with the fourth 
order system (see Figures 4.4 and 4.5). 


(4) Pressure production (term (4)), also resulting 


from density nonuniformities, primarily gives rise 
to production for u > Co and dissipetion for u < Cp: 


Since its intcsrated value over the boundary layer is 
nearly zero, this term is not thought to be of 
significance in the stability of this flow. 
To check the accuracy of the computed energy production and 
dissipation terms, each was first numerically integrated from the 


wall to the boundary layer edge and then the resultant summed. With- 


in the accuracy soupht, this sum was zero as required by 





Co 
C 


equation (2.50). Furthermore, it was found that the percent of 
positive or negative energy contribution from the production terms 
(such as the Reynolds stress) remained virtually unchanged with 
‘changes in the wall temperature level. 


It is cspecially important to note the extent of the region 


Jota 


in which Significant energy contributions take place. For the 
heated boundary leyer, in contrast to the unheated one, all 
energetics are restricted to a region very close to the wall. Out- 
Side of that region, energy terms tend to vanish, signifying that 
there is very little correlation between disturbance quantities over 


most of the boundary layer. 





CHAPTER V 


SUMMARY OF RESULTS AND CONCLUSIONS 


This investigation has examined the sixth order system of 
equations describing the hydrodynamic stability of infinitesimal 
disturbances in a heated, flat plate, water boundary layer flow. 
Throughout, comparison has been made with the fourth order 
analysis of this same problem as previously done by Wazzan, 

Okamura, and Smith [4]. 

The time-averaged (mean) ecuations of laminar flow are 
simplified by making the usual boundary laver assumptions, and by 
neglecting buoyancy forces, viscous dissipation and expansion 
energy. They are then transformed from partial to ordinary 
differential equations using Howarth~Dorodnitsyn type similarity 
variables and their resultant solutions used to provide the velocity, 
temperature, and other profiles required for the disturbance equa- 
tions. The latter are found as the linearized difference between 
the instantaneous and mean equation, after neglecting here as well 
buoyancy, dissipation, and expansion energy terms, assuming the flow 
to be locally parallel, and finally using a normal modes analysis 
to obtain ordinary, linear differential equations. The two sets of 
equations, mean and disturbance, thus derived differ from those of 
Wazzan, et.al. by the inclusion of a disturbance enerpy equation 


coupled to the other disturbance equations. This permits density, 
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viscosity, and thermal conductivity fluctuation effects to be 
considered in the stability problem. 

Both the fourth order system of Wazzan, et.al. and the sixth 
‘order system of the present investigation predict significent 
boundary layer stabilization with moderate heating, but display a 
maximum and subsequent stability decrease as the plate surface-to- 
free stream temperature difference is further increased. The 
computed stability characteristics (minimum critical Reynolds Number, 
maximum amplification rate, maximum wave number and frequency sub- 
ject to amplification, etc.) obtained from each system differ only 
slightly, however. In particular, the peaking of each characteristic 
with wall temperature (at about T EE 
neither eliminated nor significantly altered by the higher order 
system. Thus, complete stabilization apparently cannot be achieved 
by indefinitely increasing boundary layer heating. 

The insensitivity of these results to the additional con- 
sideration of thermal disturbances is attributed to the extremely 
limited region in which such disturbances. are important. Due to the 
relatively large Prandtl Number for water and homogeneous dis- 
turbance boundary conditions, the boundary layer distributions of 
thermal disturbance amplitudes "spike" very close to the wall. 

In all cases considered, the sixth order stability system is 
found to predict a slightly more stable flow than the fourth order 
system, as evidenced by larger values of the minirum critical 


Reynolds Number and reduced values of the maximum amplification rate. 
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However, it is postulated that for wall temperatures beyond the 
normal liquid range of water (with ToS 60°F), the presence of 
density fluctuations vill eventually negate the increased stability 
predicted with only viscosity fluctuations assumed, and actually 
degrade the degree of stability predicted by the fourth order 
system. That is, when therral disturbances are assumed, it is 
necessary te include all property disturbances ES obtain an accurate 
assessment of the boundary layer stability. 

The admission of thermal fluctuations also has revealed the 
existence of several energy production and dissipation terms not 
present in an energy balance equation for the fourth order systen. 
Of particular note is the disturbance shear stress production, 
ai, Zich 

dy 


the energy extracted from the mean and supplied to the disturbance 


, à source which accounts for about 10% of 


flow. In contrast, when only isothermal disturbances are assumed, 


energy is supplied completely by the Reynolds stress production, 


In summary, the fourth order system of equations used by 
Wazzan, et.ai. [4] seems to adequately describe the stability of a 
heated, water boundary layer flow in the temperature ranges 
considered. This may not be true, however, for other liquids, 
particularly those with Pr < 1, where it is conjectured that calcula- 
tions vill be much more sensitive to the influence of thermal dis- 


turbances. 
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APPENDIX A 
Detailed Derivation of Mean and Disturbance Flow Equations 
In accordance vith the discussion in section 2.1, it can 
be assumed that the flow field is describable by the set of 
conservation equations (2.1-2.3). Mean and fluctuating behavior 
are deduced from these equations in the following fashion. 
The instantaneous value of each quantity can be decomposed 


into a mean and fluctuating component. 


e = Q: (F) = Q. (x. .«) (2:56) 


where Qi is the time-averaged mean 


= Á 


Q 2 lim SR D QD. (7.1) T (A1 
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and wnere the time~averaged value of any fluctuating component 


1s defined tO De zero. 
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_ le E 
e Tu Y SG rJdr =0 


For a time average to be meaningful, mean quantitics must be 
time-independent as already denoted by the argument of 0: n 


equation (A.1). The use of time- rather than space-averaged means 
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corresponds to observed boundary layer behavior vhere disturbances 
grow or decay spatially rather than E 

Thus, substituting equations (A-1) into equations (2.1-2.3) 
‘and performing a time average, the mean equations become 


(dimensionally) 
Continuity: 
- (A.2) 


Momentum: 
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HUNE (see section 2.3), to reduce the disturbance 
differential equations from partial to ordinary ones, disturbances 
of the form Q, : 4C) exp [( 6 - e (X Co 4 Z sin 8)- cot] will be considered. 
Clearly the time’ average of this disturbance is zero, whereas a 


space averape is not (except for a neutral disturbance where eu 7 











96 


To avoid the closure problein characteristic of turbulent 
flows, it is desirable that this set of equations be solved 
independently of correlated disturbance quantities. Thus, a set of 


new nondimensional independent variables is defined 


: fi J, l 
H AN + ds = 


where E and = are length scales in the jth direction character- 
istic of mean and fluctuating flow behavior, respectively. For 
example, at a given position, L. might represent a length character- 
istic of the boundary layer development, specifically the distance 


from the body's leading edge, while 1, would correspond to the 


isturbance wave length at that point. When characteristic length 


d 
scales are comparable so that E O(1), then equation (4.5) 


reduces to Es = i Here x, must be considered a measure of both 
large and small scale variations in the jth direction, since both 
scales are in fact identical, or EE = E r= ek Such a 
situation typically aríses in this analysis where transverse 
boundary layer scales Ls and = are equivalent (exception to be 
noted later) and equal to the boundary layer or displacement 
thickness, for example, so that y = Y. 

1f it can be argued that variations in the nean motion be 


independent of small scale disturbances in each directíon (i.e., 


mean equations can be solved independently of disturhance 





Ne) 
~ 


correlations), or that the mean flow with sufficiently small 
disturbances behaves as an undisturbed flow, then Qu must be 
independent of x. (unless, of course, E = 0(1)). With this 


assumption (dimensionally) 


an Gt) = Q. (X) SS a n D (A. 6) 


and the mean equations become (nondimensionally) 


Continuity: 
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Momentum: 
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From tnese equations, several interesting observations can 
be made, 

(1) Two refcrence mean quantities are eee unspecified by 
any physical or geometric considerations: Vos 8 reference normal 
velocity; and ae a characteristic normal length scale. The 
former can be found from the continuity equation by requiring that 


a nontrivial two-dimensional solution be obtainable, so 


The latter is found by comparing the relative magnitudes of terms 
in either the momentum or energy equations. In the velocity 
boundary layer, viscous and inertial terms of the momentum equations 


are assumed equal, indicating that 





Ex = Y -y = Ke 
E 
Ly us 
Alternatively, equating conduction and convection terms of the 


energy equation in the thermal boundary layer 


Lx = Re Pr 
Ly Syste 
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Since attention herein is to be focused on the stability of a flow 
field, the length scale characterizing the velocity boundary layer 
is chosen. 

(2) For the assumption that the mean behavior be independent 
of small scale variation, Q = QG,), to be consistent, from the 


momentum equations, the criterion 
EL 
Ro me) ae 
a, Ly po 


must be satisfied, where the amplitude of all disturbances, have 
been assumed equal. If L. represents a characteristic wave 


length, then 





where eu EO >> J. Analytically, the amplitude can be assumed 
y > 


as small as is necessary to satisfy this criterion; experimentally, 
control of the amplitude is quite a diffcrent story. 
Thus, formally assuming that: 

(1) mean flow is steady; 

(2) viscous and inertial effects are equal in magnitude, so 


that attention is focused on the behavior of a velocity 
boundary layer, and 





(3) ali disturbance amplitudes are of the same order of 
magnitude and sufficiently small that correlated 
disturbance quantities need not be considered in 
solving for mean variations, or specifically 





a| «€ OL la Rey 


D 


(4) flow quantities are independent of the z coordinate 
Or the body 2S infinite in the z direction, so that 


E vz | = 
‘= s and 5 Q 


E 


©) Di. O (1) and ($) ro) but 1 /L, << 1, the 
last condition implying that the disturbance wave length 
is small compared with a length characterizing changes 
in mean quantities; 


1 se oe : 
(6) (1 fey) = O(1) and (1,/1,) z 0,95 
Mo 
Eco» 


(8) dynamic pressure is a good measure of the system pressure, 
SO ) == ; U 
S "o o0 


-1 
the mean equations to A can be Wratten as 


en A 


AA V = 0 (A.10) 
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Momentum: 

E Neu AU Y - m Op ON od 
DE. VI a = Ee N EJ | 
d sx *Y Sy) CA RR TS SY) i 
CET 2 (A.11b) 


(A.11c) 








Formally, the disturbance equations are obtained as the 

difference between the instantaneous and mean equations.  Per- 
forming such an operation with equations (2.1-2.3) and (A.3-A.5), 

the complete set of disturbance equations can be written non- 


dimensionaliy as 
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To make this set of equations tractable, the equations are 
linearized through the assumption of "infinitesimal" disturbance 
amplitudes. The linearized disturbance equations are, 

Continuity: 
B). ES NET ME ET I Io 
KS d 


to 
OU + Kr IV Ay du = O (A.16) 

















Momentum: 
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Note that the resulting equations are those of a parallel flow. 
It can be shown by this analysis that the omitted “non-parallel” 
-1/2 


terms of the boundary layer equations are of order G Re, ) : 


M y 


so that if small values of D Rey were expected, it would be 


necessary to retain these "E er, since it is expected 
that heating will increase stability [l1 - 7j, larger valves of this 
parameter are expected than for the unheated case, and the parallei 
flow assumptions will be even more accurate for this boundary layer 
flow. 

From this analysis, several interesting observations may 
be made: 

(1) To the linear order considered, there is no provision 
to account for thermal fluctuations of specific heat, second 
viscosity coefficient, or coefficient of thermal expansion (however 
slight they may be). 

(2) Additional length scales encountered in boundary layer 
stability can be easily derived from dimensional considerations. 
Consider a region very close to the boundary (viscous sublayer) with 
a characteristic velocity Ub and length E such that usb Í$ 8 


O 


and s «« n (all other conditions remaining as previously 


specified). If the boundary can be approached so closely that 


u L, ae 
e = il D e e e 
>> L = ——— , then equating viscous and inertia terms 


ob y = Ly 


ge pes 


u 
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of the momentum equation requires that 


£ Yb IT E 


of Jm 
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Incidentaliy, the pressure terms is also of the same order, 
so the full momentum equation is approximated in this region by 


Prandtl's equation [27,p 6-9] 


If another characteristic velocity, uc? and length ne are 


u 
2 
chosen so that 1 go << £ 
u x 
OC Ly 


this region, equating the inertia and viscous terms gives 
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This region corresponds to the "critical layer” where the wave 

speed and mean velocity are equal (e.s., u = Co for a two-dimensional 
flow). 

| For example, specifying Ty = 90°F, JA 60°F for a two- 
dimensional flat plate boundary layer, a solution to equations 


(A.16-A.20) with homogeneous boundary conditions is found to be 


(see Appendix E) 


Re, = 7000 
og = .122347 
Ar e D 
Ee 
Ck = .210709 
= 0 

where Cp = u at ya 5 318937. 


Thus for the length and velocity scales with the reference length 


E chosen as the displacement thickness, ô% 
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In the viscous sublayer 0 < y < ca SO that e aos > 


16.865 whereas from the continuity equation 





Ko 2r R 
Z Te WE 
UN ly = 599,57 
e E 
Se 
1 S des 
In the critical layer, y A << ya + A y so 
u 26 26 
gust 6.69/78 > =a > 3./104 and as pointed out previously, at the 
u 
critical point, Sg = 4,7459. 
e 


Thus all the length and velocity scales derived fron 
dimensional considerations seem to be in accord with derived 
solutions. Note that these two length scales, representing regions 
in which viscosity plays a dominant role, are also consistent with 
those derived from asymptotic analyses (see Mack [27, pp. 4-4,4-5] 


for example). 





Appendix E 


Thermodynanic and Transport Coefficients aud 
Their Temperature Derivatives 


B,1 Curve-Fittinp and Derivative Determination 

selection of the curve-fitting technique used in this 
investigation is predicated both on the method's capability of 
closely approximating tabulated data and for predicting high quality 
derivatives, with particular emphasis on the latter. In the ensuing 
discussion, no exhaustive analysis of all available techniques for 
achieving this objective is intended or attempted; remarks will be 


limited only to procedures considered. 


II 


Of course, a given set of data can easily be "fit" as 

closely as desired by using a least-squares polynomial of sufficient- 
ly high degree. However, if the degree is too large, small scale 
oscillations about the primary variation will be evident in the 
resulting curves due to truncation and scatter of data information. 
Utilization of these curves to generate derivatives, then, will 
naturally indicate even more conspicuous but completely fictitious 
"wiggles" in the computed derivatives. For a lower degree poly- 
nomial, not only may data be poorly approximated but variations in 
higher order derivatives might not even be seen. For both the 

method of least squares and other such empirical formulations 
considered (see relationships for p and y in section 3.2, for example) 
it is significant to note for later reference that property varia- 


tion over the entire temperature range is described by a single, 


Lu 





analytic, continuous function. Thus interpolation is not between 
"exact" tabulated experimental data points, but between approximate 
values (which presumably are within acceptable experimental error 
PRO Derivatives are then available immediately simply by 
analytically differentiating the indicated equation. 

As an alternative to least-squares polynomial curve fits, the 
method of mathenatical splines can partially satisfy the previously 
established selection criteria (assuming availability of a sufficient 
quantity and quality of data information). In general, a spline 
function is composed of piecewise polynomial arcs of degree n (where 
n is a positive integer, n » 3) which are joined at prescribed points 
such that derivatives up to and including order n-1 are everywhere 
continuous, Obviously such a function provides continuity of the 
greatest number of derivatives consistent with the use of polynomials 
of lover degree than would be required to fit all data points 
exactly by a single polynomial. Since the curvature of a mathemati- 
cal spline is most easily controlled when n = 3, or a cubic spline, 
the use of simple cubics seems to be an attractive way to interpolate 
(in the strictcst sense) experimental data points representing some 
physical relationship. As applied to this investigation, an explicit 


formula for this "piecewise cubic" is given by 


dint 
bras TE 23 A (T-T), 





DIIS 


whcre the D and the J (J > 3) distinct temperatures over the 


normal liquid range for water 


0 C = GEI Dos E x e * e e 0 o o S eq E = 100 E 


at which tabulated property data (denoted as P.) are specified 


and wnere, 


O Ton an T 


|! 


7), 
tid tor JE 


ghe t2 coefficients a, b, c, A A are then determined by 


a rare 

requiring SUE) = SR js c c cJgeaudospecitving f4rsSt or second 
derivatives et the ends of thc temperature interval, Ti or T: [20% 
The method as suggested should work perfectly if data points were 
exact. Not only would this data be exactly satisfied at the end- 
points of each polynomial serment, but the methods minimization of 
the curvature and relaxation of the overall requirement of analyti- 
city would ensure first derivatives of high quality [29]. For 
reasons intimated above, however, it is unreasonable to expect that 
data information should be satisfied exactly. Consequently, the pro- 
cedure outlined above is modified to permit departurcs from the 


prescribed data within specified tolerances while excluding unwanted 


changes in curvature. Numerically, the smoothing is accomplished by 
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ninimizing 
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where C is a positive praportionality constant and r, is the second 
J 


difference at the jth point defined by 


r = 2 J^! E E > & à [m 

EDEN x 
+1 zi i c > 
Tin ` EU je 


Note that if C = 0, then only the first sum would be minimized, and 
the procedure then corresponds to the method of least squares. For 
large C, primarily the second sum would be minimized in an attempt 
Eoreliminate undesired fluctuations [28]. Further note that in 
contrast to the first methods discussed in which temperature deriva- 
tives were obtained after smoothing data, this one computes accurate 
first derivatives while fitting the spline function to the spccified 
data. 

For the present application, the problem is thus one of 
mathematically approximating tabulated fluid property data suf- 
ficiently well that all physically realistic variation is retained, 


but spurious fluctuations are smoothed out. Note that failure 





a> 


to accomplish the latter should certainly be revealed by evaluating 
highex order derivatives. As the order of the required derivative 
is increased, the more difficult this problem is to satisfy. It is 
Guat this difficulty that determined the form of the disturbance 
equation adopted for thís analysis (see section 2.3.1). 

As an example of the methods discussed, Figure B.1 compares 
the temperature variation of the first deri DEVE of ER obtained 
from various sources by various techniques. Specifically, included 
are derivatives obtained both analytically from least-square poly- 
mental fits of Touloukian, ét.al. [30] and Kaups, et.al. [7], and 
numerically using a cubic spline fit with smoothing to the data of 
Osborne, et.al. [31]. Consistent with previous discussions, the 
latter category should indicate the more accurate general trend. 
However, it was not used for direct input to the numerical program; 
the following discussion explains some of the considerations dictat- 
ing this decision. 

(1) Even with suoothing, questionable oscillations in deri- 
vatives are indicated - specifically, those attributed to data un- 
certainty. In short, vith no unusual property-temperature variation 
appearing in the normal liquid range of water (e.g., numerous changes 
in slope and curvature), data can be well approximated by a least- 
Square polynomial of relatively low degree. 

(2) Since the amplitude of a fluctuating fluid property at 
a given nean temperature is assumed directly proportional to the 


first temperature derivative there (see equ. (2.35)), it is essential 
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l least-square polynomial 
of Kaups, et.al. {7] 
snoothed spline fit to the 
data of Osborne, et.al, [31] 
J weighted polynomial fit of 
Touloukian, ot.al.[30]! (see 
equation B,l1) 


ho 


50 100 150 200 


First temperature derivative of = E EE EE ert 


temperature as obtained using different procedures, 
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to know whether and where erratic behavier in this derivative 
occurs, Only in this vay can logical explanations be advanced for 
aberrations in eigenfunctions, Reynolds stresses, dissipation, etc. 

(3) For simplicity of evaluation and conservation of data 
storage, it is more efficient to use a low degree least-square 
polynomial or empirical curve fit. 

Similar analyses were conducted for the EOS other fluid 
properties, p, k, and y. The results suggested the property--tem- 
perature relationships for use in this stability study and are 


discussed in section B.2. 


no seeroperty-lemperature Variation for Water 

Since a major contention of this study fs that a fluid's 
property-temperature variation plays an important part in establish- 
ing its flow stability characteristics in a heated boundary layer, 
accurate values for thermodynamic and transport properties is con- 
sidered essential. The literature survey prompted by a search for 
the "best" available property information for water led to the 
following results and conclusions. 

(1) Rather than digressing on an extensive evaluation of 
data (which would entail an analysis of experimental techniques 
and the quality of information available from each, assessment of 
Systematic error Jimits in reported data, correlation of the iníor- 
mation judved most reliable, etc.), recourse is made to recent 


studies of such a nature. Althoush, as evidenced by Figure B.2, by 
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far the largest property variation with temperature is by the vis- 
cosity, equal care is exercised in specifying the behavior of each 
of the properties. In this way, discrepancies between the results 
obtained in this numerical study and those occurring in nature are 
less likely to be attributed to inaccurate property information. 
It should be noted that, of the sources encountered in the litera- 
ture survey, published data measurements found subsequent to those 
used herein agree to within experimental crror. 

(2) Since some of the "best" data available are generated by 
polynomial correlations of numerous sets of data (see below for 
E and k), concern was experienced regarding the capability of a low 
degree polynomial to be able to predict accurately the behavior of 
higher order derivatives. For example, while a second degree poly- 
nomial equation might adequately fit a set of data, the constant 
second derivative analytically derived from it may not be realistic. 
Erratic behavior in the derivetives can also be expected if the degree 
chosen is too large. Evidence of these two cases has already been 
presented in Figure B.l. Using the smoothed cubic spline as an 
approximate standard with the data set of Osborne, ct.al, [31] (judged 
to be one of the most reliable [30]), and knowing that only the 


dc, 
first derivative, : is required for this analysis, one can 


E NES 
see that the third degree polynomial of Touloukian, et.al.[30] 
adequately approximates the required behavior over the temperature 


range of interest. When second derivatives are required as well 


(such as is the case for f, p, and k), the above procedure remains 





EA 


basically the same, although more involved. For example, again 
using as an indicator the smoothed cubic spline with data recommend- 
ed by Powell [32] (also judged to be among the most reliable avail- 
‘able for k [33]), first and second derivative behavior is found to 
be best approximated when thermal conductivity is defined by a third 
degree polynomial. 

in short, when an analysis of existing data uses polynomial 
correlations which would obviously not predict behavior of higher 
derivatives needed herein, information charactcristic of that ob- 
tained from the polynomial fit is uscd instead. 

(3) Regarding the required property variation of water, the 
following information is relevant: 
Specific heat, San An analysis is made by Touloukian, et.al [30] 
of twenty sets of data for the isobaric specific heat of Jiquid 
water. For the saturated liquid, recommended values are computed 


fron 


wr 3 D 


) = 2,13974 - 9.68137x10 “T + 2.68536x10 "T 


- 2.42139x10 "1" (T in °K) (B.1) 


~} 
c (cal KR 
BS g 


BB EZ AIG N. This equation is found to fit the data 
analyzed with a mean deviation of .14% and a maximum deviation of 
1.83% (for one data set ~ all others differ by no more than .5%). 

Of particular intemest for this study, reference [30] indicates that 


c as calculated from the above cquation "should be substantially 


correct to within one percent" beJov 400^K. The ability of this 
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cubic to predict first derivative behavior has already been dis- 
cussed. 

Density po. FOU previous Works are evaluated by Gildseth, et.al: 
[34] and compared with their experimental results. Recommended 
values for the density of air free water are generated for 


0° 


IA 
I-3 


< 80°C by a modified form of the Tilton and Taylor equation 
(i.e., the inclusion of the exponential term which is negligible 


for T 240°C) 





2 | 

24. (3.9865) (T + 288,9414) OS 

nuum 1 a a nen 
(L inv C) (B2) 


to fit experimental data with a mean absolute deviation of 

ao g/ml. While the authors do specify these values to be un- 
certain in the seventh decimal] place for the lower temperatures, 
they do indicate that at least above 40°C uncertainty is in the 
sixth place and increases with temperature. Since this equation 
does reproduce the densities tabulated in the International Critical 
Tables” (35) for the range 80 s T € 100C {to within the indicated 
uncertainty), it is extrapolated in this study to cover the full 
temperature rango, 0% - 100°C. 

Dynemic viscosity, y. The experimental measurements of horosi, 

et. al. [36] are compared with ten previous, similar sets of data. 


The authors propose a correlation 
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262.37 
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to fit their data over the temperature range 20°< T <150°C with 

an average deviation of 0.17% and a maximum deviation of 0.49% at 
40°C. However, since this data seems to be characteristic of 

the mean of the other compared data for a specified temperature, an 
alternate formulation is used.  Interpolated data is peneraged 
instead by using (for the same temperature range) an equation 


recommended by R. E. Manning in reviewing the paper: 


H on SE D 
ZU C EA 1.37022(T - 20) + 8.35x10 ` = 
ee): 109 + T Loss 


to ad (20 E 1002 (cp) 
This equation is found to represent Korosi's experimental data to 
within 0.05% average deviation. 
Thermal Conductivitv, k . Of the more than sixty sets of experi- 
mental data available on the thermal conductivity of liquid water, 
Touloukian, et.al. [33] evaluate seven as being the most reliable. 
They then correlate these and generate recomnended values for 


273.16° < T < 413.16°K with the second degree polynomial 
10°: (emeu) ==" 1390.5) 4° 15.157 = 0.0190398T" (Tain ie 


This equation is found to fit the data considered with a mean 
deviation of 0.24% and a maximum deviation of 0.82%. As indicated 


in previous discussion, however, a second degree polynomial appears 


= 





122 


to be inadequate to describe the behavior of the second derivative 
in the temperature range considered. Thus, for application in this 
stability study, a third degree, least-square polynomial is fit to 
an earlier correlation of thermal conductivity data by Powell [32] 


with the following result 


k (mvatts em .K |) = - 9,901090 4- .1001982$ - 1.873892x10 ^T? 


+ 1.039570x10 q? (T in °K) (B.4) 


Values generated by this equation deviate from the preceding one by 
ne more than „23%, 

Due to the nature of the formulation of the mean and dis- 
turbance equations, it was found convenient to calculate the pro- 
duct of the density with the viscosity and thermal conductivity, 
op and pk, rather than these properties individually (sce equation 
2.20 for example). Thus temperature variations of E poo and 
p k as Vell as their temperature derivatives required for this 
analysis are shown in figures B.1 through B.5. The ratio of the 
second to the first viscosity coefficients is computed using 
equation 8.3 of reference [37, p.47] and the data of Pinkerton [14] 
and is included for reference in figure B.6 (see section 2.1). 

Since the stability characteristics are found to be sensitive 
to assumed property variation with temperature (see Chapter IV), a 
comparison is made betveen fluid properties calculated by equations 
(B.1 - B,4) (those used for this investigation) and by the Jeast- 


square polynomial fits of Kaups and Smith [7] (used by Wazzan et.al.) 
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Deviation of the latter from the former (of the form E 
Ko 
KAUPS, et.al. 
(&) IE /00 , for example, where po is the reference quantity 
an 
evaluated at |. = 32°F) is shown in Figure B.7. 
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Fipure D.2 Variations in p, EE ru, and pk with temperature, 
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Figures. 3 First and second temperature derivatives of p as a 
function of tenperature, 
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Figure B.4 First and second temperature derivatives of qu as a 
function of temperature. 
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Figure B.6 Second-to-first viscosity ratio for water as a function 
of temperature. 
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Figure B.7 Deviation between the property-temperature variation of 
Kaups, et.al [7] and that specified in Appendix B.2. 





APPENDIX C 


MEAN FLOW SOLUTION VARIABLES 


For reference purposes, the mean velocity, temperature and 
property variation through the boundary layer, as well as their 
n derivatives computed by the program are plotted in this appendix 
for the flows considered in this investigation. 

To transform this information from the n coordinate system 
to a more physical one, Figure C.1 has been included. If, for 
example, the reference normal scale is denoted as Gë Gives y= > 

O 


then 


cc CM SE 


or 


dn = yo 


ad 


In this formulation the advantages and disadvantages of using n 


(rather than y) as the independent variable may clcarly be scen. 


Advantages 
(1) Using the mean density variation, the physical step size 

is "automatically" adjusted within regions where tho temperature is 

significantly different from its free stream valuc. This means that 


although the n step size remains fixed throughout the entire in- 


tegration range, the actual distance traversed will be varied 
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constantiy within the thermal boundary layer. The larger the wali- 
to-free stream temperature difference (and so, the more severe the 
thermal gradients near the wall), then the smaller will be the 
initial steps from the wail. Admittedly, Liquid water's density 
variation with temperature is not as significant as might be desired 
for such a formulation, but for large temperature differences, it 

is still sufficient to permit sélection of a constant step size out- 
side the thermal boundary layer while ame describing varia- 
Eton within it. 

(2) If information obtained from solution of the mean flow 
equations (2.20) is to be used directly in defining the variable 
coefficients for the disturbance equations, then both must be 
expressed in terms of the same independent variable. If the latter 
is to be solved instead with y as the independent variable (to 
compensate for thc first disadvantage discussed below), then a 
constantly varying step 


x Vis 
aya EL 


| dn 
o lg e. 
1 
must be used within the thermal boundary layer. It would seem to be 
inconsistent (particularly for large temperature differences) to 
solve the mean flow equations in terms of the Howarth-Dorotnitsyn 


variables (eqn. (2.14)) and then assume the mean density to be 


constant (p = l) to solve the disturbance equations (as apparently 





Wazzan et.al. | 3-5 ] have). Such an assumption should xaise 
questions concerning the accuracy of the variable coefficients for 
the latter, especially in the wall region where density changes are 


Ie testa 


Disadvantages 
(1) Outside the thermal boundary layer, the n step size is 
LL" 


times that of the equivalent "y" one. For the temperature 


$ - er o D 
range considered mn is aivays greater than unity (E, EE E F Wren 


ó* 
D eo 2 B UC for Nox <1) so even when the displacement thickness 
is chosen as the reference length O, = ó*), a larger n step would 
actually be required than for a comparable y integration. This 
means to obtain the same accuracy, more steps would be required with 
n as the independent variable. 

(2) An extra calculation is required to transform back 
from the n to y coordinate system. 

It was decided that the advantages outweighed the disadvan- 


tages and so the problem was formulated and solved using n as the 


independent variable. 
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Fipure C.3 Second velocity derivative evaluated at the wall as 
a functlon of the wall temperature. 
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Figure C.5 Specific heat boundary layer profiles for various 
wall temperatures. 
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Fipure C6 Density bouadary layer profiles for various wall 
Fenperatutes. 
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Figure C.7 Boundary taver profiles of the first derivative of 
density for various wail temperatures, 
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Figure CO Boundary 'ayer profilcs of the density-viscosity product 
for various wall temperatures. 
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Fipure C.12 Doundary layer profíles of the density-thermal conductivity 
product for various vall temperatures, 





2| 


p 
des 
Un 


0 a: _ 
WM 
Eeer A 
60°F 
90°F 
EU 
200°F 
-.04 
-.06 N 
-.08 T = 60°F 
M = 0 
-.10 E a a E 
0 E, ‚8 > IMG 2.0 


Prourccc. l3 Boundary layer profiles of the First derivative of the 
denaity-thermal cosductlvity Pro ict. for various wall Tenperalures. 





146 


Sb 


.08 


,08 


re 


E 
O 


60 P 
90 E 


F 


150% 
TN F 


o O 9 


-,08 


T= 60°F 


Ke 


END Dope nae wee cT Ee 
0 24 des, e 126 o 


Figure C.14 Boundary layer profiles of t 
ho delay thesaad conduetivity produc 
temperatures, 


he second derivative of 
Sor various wall 





El 


147 


150 


o — 

200°F | 

150°F 

20, F 
8 50 F 
LO 
E 

T = 60°F 

a M = 0 


ER cl 


0 MÀ er | cnc TE II 
0 1 2 3 4 5 


Tiere: Boundary layer velocity profiles for various wall 
temperatures. 





148 


1.0 
T 60°F 
JD 
„8 
6 
t 
Wi 
zw 
200.7 
150 F 
; / 90°F 
l Doz 
7 
o muc o 
prp cs In Ee ee ee my [AS Ds at Sor tii Rae ny eae al mus 
0 1 2 3 A 5 
n 


fipire C.16 Boundary layer proLlles. of tie First derivative of 
velocity for various wall temperatures, 





149 


E 
To 60°F 
M= 0 
76 
» 5 
4 
M 
St 
Be $ 





0 Fa Exo sc ccc M 
0 1 2 3 4 5 
n 
piense C.I? Beundary layer profiles for the secend derivative of 


velocity for various wal) temperatures, 





150 


= 
AR DIAN DES 





o 


2 


Figure C.18 | Boundary layer temperature profiles for various 
wall temperatures. 
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Firure C,19 boundary layer prcfíles of the first derivative of 
temperature foc various wall temperatures. 
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Appendix D 


The Computer Prosram 


D.1 Program description 

The program written to solve the stability problem formulated 
in Chapter II consists of a mappines processor (BLSMAP), the main 
program (BLSTAB), and eleven subroutines, all of which are coded 
in FORTRAN V. A brief description of each, along with its role 
in the overall computational scheme and references to numerical 


procedures (where applicable), is given in Teble D.l. 
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D.2 Input 

A description of data required to successfully execute the 
program described previously is now given in Table D.2. Further, 
the required format for this data is illustrated by a sample data 
deck in Figure D.l. Note that data input is arranged so that the 
tlre. CaraelistS “only the topical control variables for the 
program; the second, all mean-flow related noe eon: and, the 
remainine, disturbance flow and stability related information. 
LE only the mean flow solution is sought, then only the first two 


cards are required; all others will be ignored by the computer. 
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Program 







d [c SQ One Ce) 


rm — d—Á— PP A 


2 1 TWALL 
2 2 TINE 

2 3 BETA 

2 4 EPSILN 
2 5 H 


l ITRMAX 


ANGLE 


Field Variable 


Program Control Parameters 


Mean Flow Paramoters 


oravit Parameters 





nn ln ee pe li 


Desert Elon 


A series of logical variables identi- 
fied in BLSTAB which controls the 
operation and output of the program. 


Wall temperature, T (specified in 


degrees Fahrenheit). 





Free stream temperature, T (specified 
in degrees Fahrenheit). 


du 
e 
Falkner-Skan parameter, M= — Io 
e 


ein 


Specified boundary layer edee para- 
meter, where KÉ =] BPSILNM. 

Step size for mean flow ecuation 
calculations, hoe It is half that used 


for the disturbance equation integra- 
tion, ha, and must be greater than or 


equal to ne /700 where ns is the 
position defined as the boundary 


layer edge. 


Maximum number of iterations per- 
mitted to compute an eigenvalue. 


An orthonormal basis is reformed if 
the angle between any two solution 
is degraded to less than this 
quantity. 


Ehe i opt m, i mn c co eri n ml A 


TABLE D.2 Description of Input Data for Program BLSTAB. 
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Program | 
1 Field p 


Card Variable escription 





Stability Parameters (cont.) 


3 S LPSEN Indicates the degree of refinement 
of eigenvalue estimate. Iteration 
is terminated when |F (ITER) 





5 


RE(EV(ITER)) — 
RE (EV(ITLR+L) ) |; and 





D 


ji . ADM (LV(ITER)) 
AIMAG (EV (ITER+1) ) | 


are all less than EPSLN. 


4 NEV Specifies which stability parameters 
are to be fixed and which are to be 
varied as the complex eigenvalue, 
EV. 

Y EV (ALFAR, OMECA) 

EV(ALFAR, ALFAT) 

EV (ALFAR, RE) 

EV(OMEGA, RE) 

EV(OMEGA, ALFAI) 

EV(RE, ALFAI) 


C^ Un E GW Pi 


1 MM Number of eigenvalues to bc computed 
(maximum of 10) 


1 REDLS (I) Reynolds number, Re^. 

2 ARDLS (I) Real part of the wave number, C 

3 ADL (T) Imaginary part of the wave number, 
Qro”. 

4 FREO(I) Frequency, 0. 


All of the four preceeding quantities 
are the best estimate of the Ith 
eisenvalue to be computed. 





| 


rnit une ee ee ër a) p ii 





Lg, fr re nn —————————————— EE nn rn nn nd nn nn 


TABLE D.2 Continued 
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D.3 Program Listing 

A complete listing of the FORTRAN V source deck used to 
‘generate the results in this investigation is given in this section. 
Note that where possible, the routines have been written in general 
so that application to systems of order greater than six is 
possible immediately with perhaps merely a change in the DIMENSICN 
statements. This is particularly obvious in the routines 
associated with boundary layer integration, ADAMS and DIFFEQ, and 
with the Gram-Schmidt Orthonormalization algorithm, GSORTH, 





BLS AR 
1 Seo Bus l Als (mE CC oS IEO- Ft (DIFF EG STABEEG-ÖRTHER-EVEST EGH FEN) ) 
oL STAB 
BL Clin Zirkel I TB I Fer), 
ee OE) pA EIS Cie ies CLO I ERE A CI O0 II ZDT ST 9,16, 351), P (13 7 3: 35 17 


1 

2 

5 Aan E EN a S RE 0(59:55] ) 

4 VISA CCO SAA 

= UMPLESN YE Vr ELUCTC+2DISTIP 
o O A AO Ak 

7 ECU CDPDEUCTOICeUJISTeFULc) 

8 eelere EECH ı CA /ELK 3 275 TU) 

g iG7(5)?:cPPCUL)/ZeGLKHSZLUPARZCZEUKSZN iHi 9 K TOT. KHAX o KMAX2 PONT Ae AR 


10 ae Ent» ARTIC (went D EUA IO TVE ESOO 

11 READ (37200) (Sale (1)+1=1+30) 

ie WRITE (67116) (GoOTECIO ^ I-15 50) 

13 IL AULA) 

i5 QUU FuKMAT (99L1) 

I5 C IFO ICA VACIA Le a OT S TRUE INE FOLEORING ACTION TS- TAKEN? 

16 COERTIOYE C2) TERMIHATE PROGRAM EXECUTION AFTE? CALCULATING AND STORING 

17 C TS EN ELO: PPRRASETERS Ati STAcIL Ely fam COE FS: 

id C uod: t2) CACCE OnE RIGC HVADUE. ege 

21 E KSE 

20 CoU DJ YRITE ETA» VELOCITY?» TEMPERATUREs AMO THEIR ETA DERIVATIVES 

21 A pP BIOPOASDEIAUGEL 

da € UE Bre Chl es TABlIPI IT ERUAFTONZCDERFTETLNTS 

a3 CNG CG? alte cuece ax EPS 

24 C TE. r ETAO ANO RATIO TOF Gle IO DISPLe« WHICKHESS 

ee Oras, sale Ice YVES 

26 C wont) wp EOVSEABMCSELEON TOOQHEANYT TXYTERSSIEPSSTZETNTLWAELTTUNESETCS 

e C cd) R TO OC eet TEPFP AMD VELOC GRAD AT WALL. AND 

e E Vial? IMCSOCSBULTISPWVARIIE AT OT TE ERROR II SOLUTE? Jee 

29 © MOTEL) ARA eee OS TAIL e Ee aA) LON Ee 

30 C Sun bc UNIT Es Fe FALSE? PAOCKAM RECORPUTES THIS IBFO 

sn C Hitt TO PROCEFLITIJO TO ETGE VALUE PROBLEM. 

GE Gere ee y eG RELATIONSHIP BETWEEN ETA AND 

a C A IS as Y A Y a E 55) 

3% E Prode.» cdi SPECIFIC I TrDcCliSITY e (CENSITVeEVISCOSITY) ¢ 

35 C Cuad CO ee dae is ETA DERIVATIVES. 

3b C Joss): e IE Le eas SIGA Z Uds Ter PTOT) IS 

Su Cure, ELO ETE Lei E e COST? Ce ilies Ae eur D FLUCTUATIONS 

A i TORNEO COMUNE ik ne el Eeler 

SY C (iert /) E ow Wi. OF Elle E EEN 

uu C lord Ye) e fev) b= Cul E m s 

4l Cau S) eee. 1 Cie EE XT WALL 
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94 
95 
90 
91 
Bd 
YY 
190 
lul 
1092 
103 
Lut 
199 
luo 
Ly? 


C 
c 
C 


(DEG 


156 


TENE PRO UE al En Sol UVECTOLSSON E DISSPOTPFERATIOM 
ES LIS: TEO EL TOS. TF FACSE USE EGUATICHNS OF WAZZANS 
OKA nid e aro EST Bam 


NC o 
S 
K nux z 201 
Berl, ul) TAE yee Bul Ar Le SI tiie 
zT pO ge laser cert) 
puc NUBE = spr ee | tet) 
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= 

Sal = Toalla 

O O 359.509 

Miia ise Wh tee rye oO 

= Pile E 

aC ieee | PE SOS. OO 
Fan 111 15579, 
Cpu TE a PUNO Ds) 
Le (MOTE (9)) ZRITE (60109) He LCETAs TWALL=459.S59rTLIIF~459,.69°PRINE? 
ll.-tErFsSlt4 

112 FORGAT t SIEP SIZE -*,D16,8/* FALWHER-SKAN PARAÁBETER zt'eF 10,67 
1 WAL TET ERATOR tFB3405.4/* FREE STARA IERMFERATURE zt',F10,.u/ 
eo! Free Stich male Woe 2'»rFY10.67*95BOUliDARY LATER EDGE DETERMINE 
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I, SO IE CI) ANO OT. NOle«l)) GO TO 1l 
o Call Meca EFS IL Pant tr | dele cE TROLS) 
Eistee EPAR KTO] 
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Ter PETE AS e A AE IE A o O E y A e LES NA A FION 
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SUSE.) 
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Ja, lei, zero 
CIO 
Ra 19989. -pudhusvcreDILSTEPADISLSETOLDHZOVAS As TALLAS TIENE 
Thee or ent lan) SIE - LEI DT TE LE TI AH I - I II I KT IK ID) 
EI pe le ar), GO 1071) 
Be ee) Tal 
Sea rt DEE ATNTELEF SAL SF TIALL I FT TEENS 
Plane wr ieee A) ELE SSA =p IMA Eb) co. tt) Go-To Y 
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EE TEE TI ORSA T TO. TS NOT THE SANC 

Vise eee ct e LE AA AA d 
STEP 
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108 
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L18 
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IT 
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jlo 
TI 
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I 
C 
ME 


Too FOS ilies 37 18EO9S 223) 
Ir ere su) ) Bist te 166 197) 
ETA Se ict fe (67108) 
Peet etal C2 soe Hite (Grito) 
lu? Foren) (4€ er Ge au FOR ECUATIOUS CF LOWELL +%+1) 
1035 Font) (7) +? FeSULTS FOR COUATIOIS OF SAZZAN?® OKAMURA» AHD SMITH 
Lo «#at/) 
Pie een (PY ees eeu Te PerRATURE FLUCTUATICUS "7 ) 
oO Elo) mae. (Seto) AMCs LTR AArEPSL™ 


ao Zr DIENEN OENB Tr et SUL VECIORS ='sF190.47 


See, ENT Drei Ton = 19/ 
et CU NERSE EE CRITERION Fun EIGENVALUES =f 2610-4) 
Hele = = 250% 
e EIER 
EE El ene USE IM EV. CALCULAETO! 
ENEE ee KE E ee RE E INe I= irn) 
249 ForkmaT (14/(4£16.0)) 
Do S Oo 
NoE Mets > PRUES 
Re == REDES (IT/ETADES 
ALFAL DIBSESCIIZCTIADES 
Q4EGA = FREOCI)ARE 
ALFAR - ARGLSCI)/ETAOL 
LSU VB PEUT 
Ie um 2) CEVA) 
If (HEV ¿Es 3) Evi1) 


CIAAPLX CALFAR OMEGA) 
CMPLX CALFAR? ALF AI) 
CMPLACALFAR + RE) 


211 11 $9 314 11 05 


IP ao. U) EVOL) CHPLXACOMEGA+RE) 
JF (NEV Ce DI EVOL) CiAPLX COMEGA? ALFAT) 
LENS NEM VEO D) EVI S Cr x CREA LEAL) 


EXE) += Mia Ue peeve) 
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Dee Oe esto CLS) (Cle O a MOTE TIO) GO TO 5 
jew ie cor) 

Pier oa. EL Ee EE EE EE EE e 
Pi lee fe sCOmRCUCTIVITY FLCOUCTUNEEOS A5SCOPEXELU) AbSULISFLU) A 
east. Felge) 

Oo + alo nPITOTY 

a SL) 

CAL rl Te 
Pos Ol Ce) ey Tak) 
Ir EE EE oe, eu Eeer? = CO Vir Do] 


Ho dn og 





Co DNO C 


FLUCTEL SR SE II + Br) File) 

PLUCTG (Crea) = G604)46P04)*Y0574)/G02) 

1: UA ALL TIGE) E. OX ERITE GEIL) Rr(FLOUCTOCOUSK) JE24 85) 
ite usto 0gx))*9523519) 

lii FORZA: (14166 16.8,/2xX+3£104,50) 

ü CORTINUE 

A e Eo) CODO. 
LO PO tenets or Eat ERPAMCE EQUATION FSR NEUTRAL DISTURBANCES. 
PLN AS) Be EECHER Eege 

[IIR Be ION DO VISCOSTTIGEBUCTUATIONSS 

Eso) POESSURE EMPRSY PRODUCTIOM 

FUrC (e) JISSIPAT IO 

FUE (Ssh) DISSIPATIO: DUE TO HOP'CODISTAHT ¡SEAN AND FLUCTUATION DEM, 
hosen) PRODUCTION OUP VO) EAN VISCOSIAY GEADIELT 

Ww«TE tier11) 

115 Porra 0/1 ENERSY PRODUCTION ANO DISSIPATIOM TERMS FOR MSUTRAL STA 
IMILITY (ALCFAI-U.0) */* KP RE STRESS PROD VISC FLUCT PROD PRES 
estribo E Eeer e E ENEE E SSC VISC GRAD PROD!/) 

Do 6 Kelek di 

L = 2uK-1 

FORCE PARTI. E) FALERRERESFREALLAY 125) 2 CO JGCLUUCUE REO $SOZECE TL 
$e COG eee (ou) IT FIRE LITER I = FUSE (19k) FC Ciel) 

Fusct2rä) AVLPAR(C3rL) 4 CREALUELUCTG(SeK) 7CONJG(T (59K) ) I ZELL IR 
Ut ALE ARE 24 AAG (CONUGCOY (19K) ) FFLUCTO CS") )) 
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SUBROUTINE MFEQN(HeEPSILNe TWALLe TINFeETADLS) 
OIMENSION ¥(15) e0Y(15) eETAEND (4) » TEST(4)»83(7» 701) »C (145 701) 


leGN(A) 
DUUBLE 
LOGICA 


GN (4) 
PRECISION Y+DY+rHr»ETA»X + Zr DELX»+DELZ»SUMODO » SUMEVEN 
L NOTE 


COMMON /BLKI/ALPHA»BETArPRINF/BLK3/G(U) »GP (4) + GPP (4) /BLK4/B ICS 


LBLKS/Ne Me KTOT eKMAX Pe KMAX2/BLKB/NOTE (30) 
EXTERNAL DIFF 
DATA ND/2/5»;(ETAEND(IJO*IZ1»2)/2.»104. / * (TEST(ID» Iz1:2)/71, » 1. E7*8/ 
l'ETEST/1. DE-12/ 
%=21,000 
Z = -1-000 
I = 1 
MARK = 1 
K= 1 

3 ETA = 0.000 
Y(1) = 0.0D0 
Y(2) = 0.000 
Y(3) =X 
Y(4) = 1.000 
Y(5) = 2 
Y(6) — 0.000 
Y(7) = 0.000 
Y(8) = 1.000 
Y(9) = 9.000 
Y(10) = 0-000 
Y(11) = 0,000 
Y(12)= 0,0D0 
Y(13) = 0.000 
Y(14) = 0.000 
Y(15) = 1,000 
CALL AUAMS(IS»HeETA«O»Yo»DYr1»DIFF) 
GO TO (4+30)» MARK 

30 B(1»K) = ETA 
B(2eK) = YC!) 
B( 30K) = Y(3) 
B(4eK) = DY(3) 
B(5;K) = Y(4) 
d(6*rK) - Y(95) 
B(7ıK) = DY(5) 

C COMPUTE MEAN N VARIATION THROUGH THE BOUNDARY LAYER 
Du 1 1=2»+ 
Gull) = di DNE VOTED 
1 Guhl) = GPP(I)*(B(6+rK)**2) + GP(I)*B(7;K) 

IF CHOTE(13)) WRITE (60101) KrGitl) + (GCI) e» GNAL ID*G(CID* GNNCI) *G(CID»I-Z 
12,3) 

101 FORMAT (I4+10E12.7) 


C COMPUTE ME 


AM FLOW COEFFICIENTS FOR DISTURBANCE EQUATIONS 


C FOR MECH PRESSURE »SET VIS2==2/3. FOR THERMO PRESSURE» SET VIS2 - 241 
VIS2 = —=2.7 56 


DVIS2 
CilsK) 
C(Z+K) 


- 0.0 
- 1.0/6(2) 
= 0G(5)/6(2) 





NOU +O 


104 


10 
IT 


12 


106 


107 


C ENGS 
13 
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C(3»K) = (2. +VIS2)*(2.*G6N(2)**x2-GnNN(2)=GN(2)*G6N(3))-=GN(2)*DVIS2 
C(4»K) 2 GN(2)*(2., *VIS2) -GNC3) 

C(5»K) = -GP(3) 

Cl(6rK) = (26t+VIS2)*(GP(2)*(35.*GN(2) -GN(3))-GPP(2)*B( 60K) )-GP(2)* 
10vIS2 

C(7*K) - -(2.*VIS2)*GP(2) 

C(6eK) - -GNCS) 

C(9»K) = -GPP(5)*B(69K)*B(39K) ~GP(3)¥*B (4K) 

C(100K) = -GP(2)*(1.+VI152) 

C(11»K) - PHINF*G(1)0/G(4) 

C(12»K) = -GNN(4) 

C(15»K) = -=2,0*G6N (4) 

C(14,K) = 24*GN(2) 


IF ((1.-Y(2))+.LE.EPSILN «AND. YC4) .LE.EPSILN) GO TO 13 


K = K+tl 

CALL ADAMSCI5seHeETA» 1» Ys» DYe A? DIFF) 

GO TO (6» 30)» MARK 

IFCETAeLTeETAEND(I)) GO TO Y 

B11 - Y(7)*x*2 + Y(9)**2 + Y(8)xx2 + Y(10)**2 

Bre = Y(7)*Y1(12) + Y(9)*Y(14) + Y(8)*Y(13) + Y(10)*Y(15) 

621 = t12 

Bee = Y(12)**2 + Y(14)**2 + Y(1L3)**2 + Y(15) **2 

Ci —OOY(2)-1.0DO)*YUCZ) + ¥CS) 4109) + Y(3)*Y(8) + YISI=*Y(10)) 

C2 = ="(1(2)-=1.0D0)+ 1112) + Y(4)€Y(18) + ¥(3)4Y(135) + Y(5)*Y(15)) 


DEM 2 Bl11*B22 =- B21*Ble 
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APPENDIX E 


The Computer Program in Operation - 


selection oBNNunsrical Parameters 
The numerical program developed to solve the stability pro- 
blem outlined in Chapters II and 111 is written in FORTRAN V lan- 
guage and evaluated on a Univac 1108 digital computer. Just as 
discussion of the complete problem has been divided into tuo 
categories (mean and disturbance), so too is the computer program: 
(1) Solution of the mean flow equations in real, 
double-precision arithmetic, and then cal- 
culation of the variable coefficients required 
for solving the disturbance equations; 
(2) Solution of the eigenvalue problem in complex, 
single-precision arithmetic and subsequent 
manipulation of the eigenfunctions to conpute 
Reynolds stresses, dissipation, etc. 
Operation of all or part of these two elements is regulated 
by a series of input logical variables. 
In addition to the parameters (or eigenvalues) which expli- 
citly characterize the forwulated stability problem (> Ors Ws 
X 
Re, and of course T? T and M), there ere also adjustable para- 
d 
meters which implicitly characterize its numerical solution (espes 
step size, orthonormalization criterion, specification of boundary 
layer edpe, etc.). Since the latter govern the ability to 
accurately and efficiently determine tne former, discussion of the 


program's operation will be developed around the optimum selection 


of these numerical parameters. 
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To Deiter illustrato ae dirtiiculties associated with 


making such a selection, numerical examples are calculated for each 


i 


parameter for the sawe typical "eigenvalue" with fixed: 


T, = 90°F, T = 60°F, M = 0, a, = 0, and Re = 7000. 


Table E.l demonstrates the sort of results obtained for the same 


= .0175, Q = 88°, n - 10, and 


d max 


fixed parameters and with: h 


E eo. 7925 (us = ,999020). 
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Comments 


j————————— TE ee 
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123/64 3212951 Including fluctuating temperature 
field but without density 


lluetudtronuss 


„121905 3.66240 Without fluctuating temperature 
field but with all mean fluid 
properties variable in dis- 
turbance equations. 


2122347 3.68283 Including all mean and fluctuating 
quantities, equations (2.37 - 2.40) 


2115509 3.360218 Without fluctuating temperature 
field but with variable mean 
viscosity only (Wazzan, et.al 
equation (1.2) d3 =~ 6J). 


. 119928 2200135 Jazzan's equation using property 
variation described by Kaups and 
Smith [7] (see Figure B.?7). 
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TABLE E.1 The Influence of Included Disturbance Quantities and 
Chosen Property Variation on an Eigenvalue Calculation. 
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E.l Boundary Layer Edge, n 


and 
max do 


UA A A es m nn - 


Two distinct length scales must be specified to accurately 
Solve the eigenvalue problem: ee the finite value of the 
independent variable, n, at which the mean flow's asymptotic 
boundary conditions (eqns. (2.21la)) are assumed to be approximately 


satisfied; and Mg (n, « n ), arbitrarily selected boundary layer 


max 
edge where asymptotic disturbance flow conditions are satisfied. 
For ccmputational purposes, there is no reason to sugrsest that the 
two need be the same. While the former is associated with obtain- 
ing accurate mean velocity and temperature profiles, the latter 
determines where these profiles may be truncated for the eigenvalue 
problem without appreciably degrading the accuracy of its solution. 
The practical capability for determining the proper Ne 
however, is contingent upon utilization of an efficient means for 
calculating the initially unspecified, mean-flow, wall boundary 
conditions, u'(0) = F"(0) and H!'(0). The effectiveness of the 
particular method used in this investigation [18] for finding these 


conditions is demonstrated by the example in Table E.Z2 where 


T = 90°F, tT = 60°F, M= 0, h_ = .02, n._ = 10. 
W co m nmax 





3.59 


em A 


Iteration ng pom BON 
Guess 120 -1.0 
il 2 0. 63890603 -0.80589459 
2 10 049 Tomos -0.71569664 
3 10 0.46637891 -0.71328643 
| 4 10 0.46639972 ASA 
5 10 0:460939971 +0. 71321361 
6 10 0.46639971 0.7182 bool 
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TABLE E.2 Convergence History of Unknown Mean Flow Boundary 
Conditions, © (0) and h (0S 


With 0.55 seconds required for the first iteration, approximately 
2.78 seconds/iteration thereafter, and six iterations, both con- 
vergence to values of F'(0) and H'(0) with better than eight-place 
accuracy (for the specified fluid property variation, of course), and 
calculation and allocation to storage of all mean-flow coefficients 
needed for solution of the disturbance equation can be accomplished 
in less than twenty seconds (see Table E.3). 

Final selection of De is determined by observing the value 
for which an increase in this parameter fails to make a discernable 
change in F'"(0) and H'(0). It is reasoned that values of N 
greater than this one will not alter the mcan flow profilcs (thus, 
the variable coefficients for the disturbance equations) and so 


must represent an unnecessary expenditure of computational effort 


and tive. 





190 


n 
a 


The results of such an Wees Search for the mcau wall 
temperature and velocity gradients considered in this study are 
presented in Table E.3. In all cases, F'"(0) and H'(0) were un- 
changed for Dam 10. Thus, noting that the time required to 
solve the mean flow equations and compute disturbance equation 
coefficients varies only slightly over the Tmax range considered, 


a 


accuracy rather than econony governed the selection: 


This value ensures satisfaction of asymptotic conditions (eqn. 

(8 4) sto order a0 19), Due to its importance in calculating the 
eigenvalue, the influence of Hit on Dos is also displayed in 
Table E.3. As a basis of comparison, it may be noted that the 
value of F"(0) for the unheated, flat plate case agrees well with 
Howauen sevalue ol FO) = — 33206 T9, p: 129]: 


Tbe primary considerations for selecting were founded in 


"e 
attainment of an accurate eipensolution and reduction of running 
time and storage requirements (known a priori to be relatively 

large when eigenfunction recovery is required [20]). Although the 
time restriction could be relaxed somewhat if the information were 
needed badly enough. storage limitations (130,000 for the Univac 
1108 computer) could not. To evaluate the effect on the eigenvalues 


of using a "truncated" mean flow solution, the example shown in 


hap le e. S vas Calculated vith: Ue = 90°F, 











F"(0)-u (0) u' (0) ng” E (defined Time 
by Eqn. Required 
(sec) 
Ty = 60°F 
6 23255559 iets GMS o sica 
| 7 | .33202439 1.7203045 .43557x10 1500072 
8  .33205907 1.7207617 .13211x10 ^ 16.7630 
9 — .33205738 1.7207868 .14727x10 l^ 17.7298 
10  .33205735 1.7207876 .60392x10 |? 18.8970 
ll  .33205733 1.7207876 .91101x10 ^ 20.9158 
T = 90°F 

6 —.46680267  -.71341605 1.5112973 .19198x10 ? 14.6256 
?  .46642702  -.71322750 1.5142917 .12898x10 ^ 16.1574 
8  .46640087  -.71321419 1.5145365 .31727x10 ll 15.2932 
2807565997209 71391562) 1.5145486 286951072 17.0484 
10  .46639971  -.71321361 1.5145491 .95693x10 ?’ 18.9308 
E 0097 1929960081 0125491 999252193102 5299195 2002 

„= 150% 
6  .74345286  -.78390563 1.2098:09 .35870x10 ? 14,5934 
7  .74320905  -.78382177 1.2111477 .17598x10 ? 14.5190 
Bin 7o210404 78381462 162112361 .31693%10 15.8284 
9  .74219379  -.78381642 1.2112399 .21000x10 |^ 16.9238 
10  .74219378  -.78581642 1.2112400 .68194x10 |? 17.7392 
11 nenn  -.78281642 1.2112400 .17341x10 |? 19.3116 
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TABLE E.3 Influence of n, . on computed values of F'"(0), H'(0), 


and ne^. 











T 
w 


6 
7 
8 
2 
10 
11 


T = 60°F = 0. hi 
= TY), 


= 200°F 








. 96336141 
263188592 
290317729 
223051 7605 
290517604 
~ 96317694 


TABLE E.3 “Continued 


branch without density 


BUS 
7.00 
6.70 
5. £0 
4.68 
3.14 





EE ————À— — M — "9 O 


M-0,h 


"e 
.999959 ` 
.999903 
. 999038 
. 990055 


980150 





197 


nn fe mm oe e 





a ee ee eee HA 








-.83923358 1.0362211 .11659x10 9 12.3970 
-.83918394 1.0359269 .47597x10 ! 13.5738 
-.83918114 1.0369723 .71388x10 l^ 14.8100 
-.83918104 1.0369741 .39403x10 7° 16.5640 
1559915104. 85055975999 5085/2510. 22 18.070 
-.83918104 1.0369741 .61905x10 1° 18.5082 | 
- .02 
hos 02, 0 = 455, Re¿* = 7000, a, = 0 (lower 
fluctuations). 
Soto AAA 
f Steps normali- Tíme/ Go Se 
_..-__._._ 2ations Iteration 0 ^ » 0 
351 102 li.i514 .123762 3.72956 
Bon 102 O EE 
291 99 9.5320  .123817  3.72989 
235 96 1565708 291245669099 :579520 
158 91 5.1472 .134641 3.83808 


TABLE E.4 Influence of the Arbitrarily Specified Boundary Layer 


Edge, Nez ON the eigenvalue computation. 
3 


Note that consideration of additional segments of the 


boundary layer does not appreciably change the number of ortho- 


nornalizations required, suggesting that most of them are required 


in the region near the wall. 


It does, however, require more 





integration steps and so more time per iteration. Unfortunately, 
due to storage and time limitations, it is more difficult to as- 
-cribe a definite value to Ns (also the case for the disturbance 


equation step size, ha) than for i . Even using the maximum 


i 
e 


Clas 


allocated storage (ha = .02 and a 7.0), varıalaons in Te 
eigenvalue axe still slightly greater than their indicated six 
place A 

Since even four place accuracy is adequate to clearly 
delineate between calculations with and without temperature 
fluctuations (see example in section 3.4), it is sufficient to 


specify 
e 


as the position at which the boundary layer may be truncated for 
stability calculations. 
Greater accuracy than that indicated in Table E.4 is 


restricted by both storage and the interdependence of Ns and ha TU 


the program; when h, must be reduced, so too must Ns according to 


d 


the relationship h, > ng / 350. 


d 





E.2 Step Size, h_ 

Perhaps no other parameter more critically governs the 
accuracy and efficiency with which a solution is numerically cal- 
culated than does the step size. If it were chosen too large, the 
fourth order Runge-Kutta and Adams~Houlton integration schemes 

t 
(with truncation errors proportional to h^") would be basically 
inaccurate, rapid changes of the dependent variables (such as in the 
thermal boundery iayer for large Pr) would be obscured, and (most 
importantly for this analysis) between consecutive steps, for a 
given oRe, linear independence of the solution vectors may be 
destroyed beyond accurate restoration by the orthonormalization 
procedure, independent of the $ specified (see section 3.4.3). If 
chosen too small, not only would integration time and storage be 
prohibitively large, but error propagation, enhanced by the in- 
creased number of mathematical operations required, would contami- 
nate the solution sooner and thus require more frequent ortho- 
normalization within a given interval (see section 3.3), 

Thus, to solve the disturbance equations, the selected h 
must be very closely related to the specified values of ake and Q. 
Since accurate convergence is ensured only by giving the ortho- 
normalization procedure a chance to work effectively, as aRe 
increases, h should be decreased and 9 increased. In mapping out 
aestabrility contour, the point at which convergence difficulty can 
be anticipated js preceded by a limited ranse of smaller aRe in 


which the number of iterations required for convergence is greatly 





increased ana more importantly, the minimum angle encountered be- 
tween solution vectors drastically drops to less than 1°. At or 
beycnd this point, convergence could be obtained only be decreasing 
the step size. It is interesting to note that similar difficulties 
were encountered by Mack [24, p.283] (in a scheme not using 
orthonormalization), but at much smaller values of oRe. 

In contrast, for the numerical values calculated using the 
fourth order system of Wazzan, et.al [3]. no such difficulty was 
found, One possible explanation can be advanced by considering 
the examples given at the besinning of section 3.4. For the sixth 
order system (that is, the problem as formulated in this study but 
without density fluctuations), the minimum angle encountered during 
the boundary layer integration was 36.40%, whereas for the fourth 
order system (Wazzan's equations) it was 78.21". Since the region 
ii which significant thermal fluctuations exist is much smaller 
than that of the velocity fluctuations, it can be expected that the 
step size accordingly must be smaller for the sixth order system 
to maintain the same degree of orthogonality as for the fourth. 

Since the Runge-Kutta integration scheme requires infor- 
mation at half-step intervals, the mean equations are integrated 
with a step size, n.o half that of the one used to solve the 
disturbance equations, D (PUES D. - IR This permits the 
program to calculate the required variable, mean-flow-dependent 


eoetilicients for the latter exactly for all points. 





F2 
M 
ON 


An example illustrating some of the points discussed above 


ts piven in tables E.5 and E£.6 


n = 10, u, =. a5 


| max Â 


without density fluctuations). 


u (9) 
46639971 
46639971 
‚46639972 
46639973 
46639977 
46639982 





















for: 10290 7, T. = 
W E 
Res” = 7000, Ay 
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H (0) 
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2309501 
219321301 
2021521559 
E0927 
cO ESSE 
EDS? 


ele 
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gon) 
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1.5145489 
ISS 


60°F, M = 0, 


: 0 ( lower branch 









TABLE E.5 Influence of Step Size, hos on the Unspecified Wall 


Temperature and Velocity Gradients, u'(0) and H (0). 


m i i i ll o A P i A APRI p t M 


Time/ 

ha Steps Iteration( 
701673 347 11.1042 
. 02 291 9.5320 
203 194 6.4768 
04 146 5,0016 
203 117 4.0280 
.06 98 39590 


Oo Å- á= rm a 


TABLE E.6 Influence of the Ste 


Eirenvalue. 


#Orthonor= Go 
Sec) malizations 
103 2129910 
Do 2129917 
94 2123082 
78 2129797 
67 2975 
59 223 





D oiz; hg’ on the 


Computed 


"ST 

3.729390 
Oy 12909 
SU aT 
gon y 
3.72834 
Do UD 





From these calculations (oke = 866.6), it appears that at least 
four place eigenvalve accuracy is obtainable by using A < 409% 
Note that this choice also ensures eight place accuracy of the 
mean flow solution (for the specified property variation). For 
results given in Chapter IV, a step size of ha = .02 or smaller was 


used. 


E.3 Orthonormalization Angle, 9. 
As previously indicated (section 3.3), the orthonormalization 
criterion used for this study is a specified mininun angle per- 
mitted between complex solution vectors. (A comparison of relative 
magnitudes of the vectors also could have been used [39]). If 
this angle were chosen too large (e.g., 90? as for Wazzan, et.al. 
[3]), aside from an increased integration time, the extensive 
matrix M required fox the orthcnormalization procedure 
could be expected to introduce errors in the solution [39 ]. 
Alternatively, if it were chosen too smali, then linear independence 
of solution vectors could be so degraded that insufficient resolu- 
tion would remain to accurately reestablish an orthonormal basis. 
Since the parametric froup aRe appears exponentially in 
the exact solution to the disturbance equations outside the boundary 
layer (see equations (3.11 - 3.15)) or inversely as a coefficient 
multiplying the hieest order derivative term (see equations 
(1.2 - 1.2)), it can be anticipated that as wRe increases, so too 


should the rate at which the solution vectors linear independence 





mE 


| 


P a <a 


A 


p 


degenerates. It follows, then, first that for the specified Q, 
the number of orthonormalizations required must also increase. 
Secondly, as 2 is increasing to accommodate the larger value of 
aRe, so too should the percentage of orthonormalized intepration 
steps, ranging from zero when 2 = 0° to one hundred when 2 = 90°, 
These observations are substantiated by Figure E.l and Table E.7. 
Care must be exercised to ensure that Q is chosen suffi- 
ciently large that the resolution of solution vectors is not lost 
betveen successive integration steps just before the orthonox- 
malization criterion is violated. That is, if the smallest angle 
between solution vectors at one step is greater than Q, but less 
than 2 at the next, this minimum angle must still be large enough 
to enable reformation of the orthonormal basis. When degeneration 
of orthogonality is sufficiently rapid to require orthonormaliza-— 
tion at consecutive steps, it can be expected that the minimum 
angle between solution vectors (and thus convergence capability) 
would cease to be a function of 2 and become one of the step size 


instead, as indicated in Figures E.1 and E.2. 
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9.0 
T = 909r 
8.0 T, = 60 F 
M = 0 
Ds ,999 
Q 45? 
2m 
1.0 Re _ n lower branch 
SE r=0 
0.0 
90 
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bipure E.2 Influence of the step size on the minimum angle 
encountered between solution vectors on the last eipenvalue 
Ptercttion for a specified, 









Number of 





Tre 

5 2179330 p DE 44 1.3507 
10 2079350 3. 70952 50 105990 
20 -179330 2. 10852 74 1.6808 
30 O 3,70852 Z9 1.6808 
40 21710530 5-705592 83 1.68608 
59 211293930 USO Se ` 1.6608 
60 “179330 270992 93 1.6808 
70 1179330 3210852 100 1562. 7:33 
60 29530 3.70852 H? ee 
88 272330 3: 7/09533 140 12.3308 
90 700g DS 146 not computed 


: 4 
where Re = 10 , «a 


SEE ccr. u, 


S = .999 (upper branch, r # 0) 


i 


TADLE E.7 Influence of Orthonormalization Angie, 2, on Computed 
Eigenvalue. 


On the basis of these calculations, it would seen that al- 
though the minimum angle between solution vectors is profoundly 
affected by the specified orthonormalization angle in some 
regions, the computed eigenvalue is not. Further, in the 
intermediate range of &, there is a region in which the number 
of reauired ortnonormalizations is significantly less than for 
Q c 90° and where this number varies only slightly with Q. 
Naturally, such a reduction translates as a significant savings 
in required program running tire. For example, for the Re = /000 


(upper branch) calculation shoun inu Figure E.1, the time reauired 





202 


to complete eigenvalue convergence (in seven iterations) is re- 
duced from 36.9 to 34.4 seconds by specifying 2 equal 50° rather 
than 88°, still with almost six-place accuracy. Thus, selection 
of 2 apparently can be based primarily on computational econony; 
that is, it is chosen to reduce the number of computer manipulations 
nd progran storage requirements. 
Since these results are typical of others found during tne 


course of this study, it is fclt that the selection 


De au, 


represents a good balance between accuracy and efficiency. To note 
in passing, this choice is also compatible with the results of 


Gersting and Jankovski [20]. 


E.4 FEirenvalue Estimates, a, u, and Re 

As indicated in sectíon 3.23, three eigenvalue estimates are 
required to start either of the root-findiüp schemes used in this 
Study (ice: bhariof Wazzan, et.al: jor Muller). To trace out a 
contour in parameter space, the first Ue estimate, evi, is 
taken from the lower branch of the stability curves of Wazzan, et.al. 
[4] near the minimum critical Reynolds number (see for example 
Figure 4.4 ). This procedure is employed for two reasons. 
Tirst, since the present investigation represents a departure Írom 


the Wazzan results due to the inclusion of a fluctuating, temperature 


field and consideration of all fluid property variation, it would 








seem that these should give the best available initial estimate. 
second, as the difficulty of solution increases vith increasing 

cRe (due to problems associated with degeneracy of the orthogonality 
of solution vectors, etc.), an eigenvalue on the new contour is 
chosen for which solution should be easiest (see also discussion by 
Mack [24, p. 280]). 

The second and third eigenvalue estimates are arbitrarily 
chosen as ev) = erg evi and ev, = T5005 Vir oque y Car: Dreh 
the three initial eigenvalue estimates thus specified, the three 
corresponding boundary layer admittances, Yo can be calculated. 

The plane fitting technique of Vazzan, et.al. is used first, 
conditionally followed by tbat of Muller. That is to say, only 


when all the following criteria are violated will the latter be 


utilized: 
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(2). Pirom the denominacor of a 
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toa < 10 
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(3) Muller's method was not used in computing the 
preceding estimate. (It was found that consecutive 
application of this method gave poor convergence.) 


In general, this scheme worked well, ensuring six place 
accuracy within five to cight iterations (depending on the 
accuracy of the initial guess). At first, to ensure uniform con- 
vergence (after the first three arbitrary estimates, ev, > eva), 
an estimated "eigenvalue" obtained by either metlod was first 
averaged with the most accurate one previously calculated (i.e., 
that one which gave the smallest value of IY D before reevaluating 
the stability equations. However, since this procedure was found 
to increase (by as much as 50%) the number of iterations required 
for convergence, it was modified to average only when the condition 
of uniform convergence was violated (i.e., when the condition 
Io 1) > Yo! was not satisfied). When Muller's method was used 


alternately as discussed, it was found to be very effective, usually 





refining the "cigenvalue” cstimate so as to reduce za by several 


orders of magnitude., Thus, the net result is an iteration scheme 
which ultimatcly but sporatically converges, often in a nonmontonic 
fashion. 

An example of the convergence is given below for DE 90°F, 


T = 60°F, M= 0, 2 = 50°, h, = .04, u, = .999, Re = 7000, o 


6 = Us 


I 


-6 
Na ag O o AO "(Lower branch including density 
R ? Te 


Miao uations r y 0): 
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lteration 





















Cuess 1 .081212293 .017239455 .432x10 7l 
Guess 2 082024415 017411849 .342x10 5 
Guess 3 . 081618354 017584244 .139 
1 081367993 .017209441 .176x10”1 
2 .080728449 .017008234 .332x10 ~ 
3 080769324 017018352 323x107 ^ 
4 .080769138 .017018221 Eon 
next estimate 080769140 2017015221 





TABLE E.S Convergence History of an Eigenvalue Search 


In keeping with previous discussion of obtainable and re- 


quired accuracy, iteration was terminated when both 





IY | 
O 
and 
ev, ev, 
1 - A and l -— - < 10 > 
TLR (j+1)1 


This criterion usually ensured "eigenvalue" accuracy to six places 


and "zero" wall boundary conditions of order (107^ OL less. 
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